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Notation 


A,  B  State  and  control  matrices 

A(t)  Open  loop  stability  matrix 

Aq  Stability  matrix  (constant  approximation) 

A,  B  A  and  B  matrices  in  Fourier  reformulation 
A(t)  Closed  loop  stability  matrix 

B  Number  of  blocks  in  the  Harmonic  Transfer  Function 

B(t)  Input  (control)  matrix 

Bq  Input  (control)  matrix  (constant  approximation) 

C,  D  Output  and  feedthrough  matrix 
C(t)  Output  (measurement)  matrix 

Co  Output  (measurement)  matrix  (constant  approximation) 

Ct  Rotor  thrust  coefficient 

C.  T>  C  and  D  matrices  in  Fourier  reformulation 

D(s )  Denominator  of  the  system  transfer  function 

D(t)  Direct  input /output  matrix 

D0  Direct  Input/Output  matrix  (constant  approximation) 

Ek  Feedthrough  matrix  of  time-lifted  reformulation 

Fk  State  matrix  of  time-lifted  reformulation 

G(a,  k )  Input-output  transfer  function 
Gk  Control  matrix  of  time-lifted  reformulation 


G(s)  System  transfer  function 

H(s)  IBC  compensator’s  transfer  function 

Hk  Output  matrix  of  time-lifted  reformulation 

Hi(z,  s )  Sampled  transfer  function 
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m  Number  of  control  inputs 

M  Monodromy  matrix 

k  Sample  number,  k  =  dT  +  s,  and  discrete-time  variable  associated  with 

T  (1/rev  dynamics) 

K  Gain  of  HHC  controller 

t  Number  of  periods  required  to  reach  k 

Mt(k)  Periodic  Markov  coefficients 

N  Number  of  rotor  blades 

n  Number  of  measured  outputs 

n  System  order 

ns  System  order 

ns  Number  of  acceleration  (output)  samples  per  revolution,  ns  =  T / P 

N(s )  Numerator  of  the  system  transfer  function 
p  Number  of  measured  outputs 

P  Acceleration  (output)  sampling  interval,  P  =  T /ns 

r  Relative  degree 

r  Tuning  parameter  in  HHC  performance  index 

R  Zero  dynamics 

R  Control  weighting  matrix  in  HHC  performance  index 

s  Sample  number  within  a  period 

T  Number  of  samples  in  one  period,  also,  period  of  one  rotor  revolution 

T(j  HHC  gain  matrix 

t  Time 

u  Control  (or  input)  vector 

Pilot  or  flight  control  system  input 
vlhhc  HHC  input 


v(fc),V(cr)  Generic  signal  and  its  ^-transform 


W(a) 

W 

x 

y 

y nc  y m 

Z 

P 

r 

v 

Oo 

01ci  &ls 
0kc y  0ks 

n 

<p 

$ 

tp 

^(0) 

a 

a 

uli 
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Heave  velocity  and  acceleration  of  the  CG 
Frequency- lifted  transfer  function 
Output  weighting  matrix  in  HHC  performance  index 
State  vector 

Output  (measurement)  vector 

iV/rev  cosine  and  sine  harmonic  of  i-th  acceleration  measurement  (output) 
Time  shift  operator  (by  T  samples) 

Blade  flapping  angle 
Instantaneous  interaction  matrix 

Discrete-time  variable  associated  with  P  (acceleration  sampling) 

Collective  pitch 

Lateral  and  longitudinal  cyclic  pitch 

k/rev  cosine  and  sine  harmonic  of  pitch  control  input 

Advance  ratio 

Complex  number,  cf)  =  exp(j2n/T) 

Transition  matrix 
Transition  matrix  from  r  to  t 
Azimuth  angle  of  reference  blade,  ip  —  Of 
Transition  matrix  at  the  end  of  one  period 
Rotor  solidity 

Time  shift  operator  (by  one  sample) 

Fundamental  rotor  lag  frequency 
Rotor  angular  velocity 


Subscripts  and  superscripts 

(  )c  Quantity  in  the  controller  model 
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(  )  f  Quantity  in  the  harmonic  analyzer  model 

(  )  h  Quantity  in  the  helicopter  dynamic  model 
(  )uft  Time-lifted  quantity 

(  )z  Quantity  in  the  zero-order- hold  model 
(~)  Discrete-time  counterpart  of  continuous- 

time  signal 
Abbreviations 

CG  Center  of  mass  of  the  helicopter 

EMP  Exponentially  Modulated  Periodic 

HHC  Higher  Harmonic  Control 

HTF  Harmonic  Transfer  Function 

IBC  Individual  Blade  Control 

LTI  Linear  Time  Invariant  (or  constant-coefficient)  model 

LTP  Linear  Time  Periodic  (or  periodic-coefficient)  model 

MCT  Multiblade  Coordinate  Transformation 

MIMO  Multi-Input  Multi-Output 

NMP  Nonminimum  phase 

PHHC  Periodic  Higher  Harmonic  Control 

RHP  Right-Half-Plane 

SISO  Single-Input  Single-Output 

/rev  Per  rotor  revolution  (frequency) 
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Chapter  1 


Summary 


This  Final  Report  begins  with  a  summary  of  key  accomplishments  for  the  entire  project. 
They  include:  (i)  the  first  systematic  study  of  the  effect  of  zeros  in  rotorcraft  aeromechanics, 
(ii)  the  first  study  of  the  effects  of  closed-loop  HHC/IBC  on  the  aeroelastic  stability  of  a 
helicopter  rotor,  (iii)  a  presentation  to  the  rotorcraft  community  of  techniques  for  the  analysis 
of  systems  with  periodic  coefficients  developed  in  other  areas  of  engineering,  and  (iv)  the 
first  study  of  the  effects  of  closed- loop  HHC/IBC  on  the  aeroelastic  stability  of  a  helicopter, 
considered  as  a  discrete  system.  Each  of  the  accomplishments  and  their  significance  are 
briefly  discussed. 

The  next  chapter,  Chapter  3,  describes  the  role  of  the  zeros  of  input-output  models 
for  coupled  rotor- fuselage  systems,  modeled  as  linear  systems  with  constant  or  periodic 
coefficients.  The  zeros  do  not  affect  open  loop  stability,  but  they  can  affect  closed  loop 
stability,  and  may  trigger  closed-loop  instabilities  with  high  feedback  gains.  The  constant 
coefficient  approximations  obtained  through  a  multiblade  coordinate  transformation  usually 
give  accurate  pole  predictions  for  advance  ratios  of  up  to  n  —  0.3  —  0.35,  but  the  zeros  are 
reasonably  accurate  only  up  to  about  /i  =  0.2.  The  calculation  of  the  zeros  can  have  much 
worse  numerical  properties  than  that  of  the  poles  because  of  large  ratios  between  largest  and 
smallest  zeros.  Therefore  the  techniques  normally  used  in  Floquet  stability  analysis  often 
prove  inadequate.  The  study  also  shows  that  the  lag  poles  are  not  canceled  by  corresponding 
zeros,  therefore  the  heave  response  to  collective  is  affected  by  rotor  lag  dynamics,  especially 
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at  higher  thrust,  and  for  stiff-in-plane  configurations.  This  transfer  function  has  several  zeros 
close  to  the  imaginary  axis  or  with  positive  real  parts.  Their  number  and  location  depend  on 
the  configuration,  but  at  least  one  nonminimum  phase  zeros  is  present,  and  at  every  speed. 
This  limits  the  benefits  achievable  from  active  control  systems. 

The  following  chapter,  Chapter  4,  describes  the  development  of  a  state  space  formulation 
for  a  Multi-Input  Multi-Output  (MIMO)  Higher  Harmonic  Control  (HHC)  system.  The 
development  starts  with  a  simple  state  space  derivation  for  the  continuous  time  form  of  a 
Single-Input  Single-Output  (SISO)  HHC  compensator  with  input  and  output  at  the  same 
rotor  harmonic;  then  the  same  approach  is  extended  to  the  case  of  different  harmonics  in 
input  and  output,  which  result  in  a  periodic  SISO  HHC  compensator;  finally,  that  result  is 
generalized  for  the  derivation  of  the  state  space  form  for  a  MIMO  HHC  controller.  The  chap¬ 
ter  contains  results  of  a  numerical  investigation  into  the  performance  and  stability  properties 
of  a  closed  loop  HHC  system,  implemented  in  the  rotating  system,  based  on  a  simulation 
study  of  the  coupled  rotor-fuselage  dynamics  of  a  four  bladed  hingeless  rotor  helicopter. 

The  results  show  that  the  HHC  controller  is  very  effective  in  reducing  the  4/rev  acceler¬ 
ations  at  the  center  of  gravity.  The  percentage  reductions  obtained  in  the  simulations  are  in 
excess  of  80-90%.  The  vibration  attenuation  occurs  within  5-7  seconds  after  the  HHC  system 
is  turned  on.  This  is  equivalent  to  a  frequency  of  around  1  rad/sec,  which  is  a  frequency  at 
which  flight  control  systems  and  human  pilots  tend  to  operate.  Therefore,  the  interactions 
and  potential  adverse  effects  on  the  stability  and  control  characteristics  of  the  helicopter 
should  be  explored.  The  HHC  problem  is  intrinsically  time-periodic  if  the  HHC  inputs  in¬ 
clude  frequencies  other  than  the  frequency  one  wishes  to  attenuate.  This  is  true  even  if  the 
rest  of  the  model  is  assumed  to  be  time- invariant.  In  these  cases,  the  closed-loop  stability 
results  obtained  using  a  constant  coefficient  approximations  may  be  incorrect  even  at  lower 
values  of  the  advance  ratio  /i,  where  constant  coefficient  approximations  of  the  open-loop 
dynamics  are  accurate. 
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Chapter  5  describes  in  great  detail  some  basic  concepts  in  the  treatment  of  discrete 
systems  with  periodic  coefficients.  In  particular,  the  chapter  describes  the  derivation  of 
transfer  functions  for  such  systems,  which  are  the  foundation  for  the  methods  to  convert 
systems  with  periodic  coefficients  into  equivalent  systems  with  constant  coefficients. 

The  following  chapter,  Chapter  6,  describes  how  to  convert  the  periodic  coefficient  rep¬ 
resentation  of  a  discrete  linear  system  into  a  constant  coefficient  model  that  has  the  same 
stability  characteristics  and  is  equivalent  from  the  input-output  point  of  view.  This  is  the 
key  ingredient  for  the  stability  analysis  of  a  helicopter  with  a  closed-loop  higher  harmonic 
control  (HHC)  system.  In  fact,  typical  HHC  implementation  contain  discrete  elements  and 
multiple  sample/update  rates.  The  extraction  of  a  linearized  discrete  model  from  the  con¬ 
tinuous  one  is  described.  The  Floquet  characteristic  multipliers  and  characteristic  exponents 
are  the  same  for  the  discrete  and  continuous  systems,  and  therefore  the  open  loop  stabil¬ 
ity  characteristics  of  the  two  systems  are  exactly  the  same.  Time-lifting  converts  a  system 
with  periodic  coefficients  into  a  larger  system,  but  with  constant  coefficients.  Time-lifting  is 
based  on  two  properties  of  discrete  periodic  systems,  namely:  (i)  the  response  of  the  periodic 
system  can  be  modeled  as  the  combination  of  the  responses  of  multiple  constant-coefficient 
systems,  one  per  sample,  and  (ii)  using  shifting  operators,  it  is  possible  to  “move”  all  the 
system  inputs  at  the  same  sample  in  multiple  periods  to  that  sample  in  a  single  period. 
Other  discrete,  constant  coefficient  reformulations  are  available,  both  in  the  time-  and  in 
the  frequency-domain.  These  reformulations  can  be  useful  for  theoretical  studies  and  prac¬ 
tical  applications.  Numerical  examples  for  a  simple  helicopter  blade  model  are  provided  for 
transfer  functions  and  time-lifted  reformulations. 

The  last  chapter,  Chapter  7,  presents  an  aeromechanical  closed  loop  stability  and  response 
analysis  of  a  hingeless  rotor  helicopter  with  a  Higher  Harmonic  Control  (HHC)  system  for 
vibration  reduction.  The  analysis  includes  the  rigid  body  dynamics  of  the  helicopter  and 
blade  flexibility.  The  gain  matrix  is  assumed  to  be  fixed  and  computed  off-line.  The  discrete 
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elements  of  the  HHC  control  loop  are  rigorously  modeled,  including  the  presence  of  two 
different  time  scales  in  the  loop.  By  also  formulating  the  coupled  rotor-fuselage  dynamics  in 
discrete  form,  the  entire  coupled  helicopter-HHC  system  could  be  rigorously  modeled  as  a 
discrete  system.  The  effect  of  the  periodicity  of  the  equations  of  motion  is  rigorously  taken 
into  account  by  converting  the  system  into  an  equivalent  system  with  constant  coefficients 
and  identical  stability  properties  using  a  time  lifting  technique. 

The  most  important  conclusion  of  the  study  is  that  the  discrete  elements  in  the  HHC  loop 
must  be  modeled  in  any  HHC  analysis.  Not  doing  so  is  unconservative.  For  the  helicopter 
configuration  and  HHC  structure  used  in  this  study,  an  approximate  continuous  modeling 
of  the  HHC  system  indicates  that  the  closed  loop,  coupled  helicopter-HHC  system  is  always 
stable,  whereas  the  more  rigorous  discrete  analysis  shows  that  closed  loop  instabilities  can 
occur.  The  HHC  gains  must  be  reduced  to  account  for  the  loss  of  gain  margin  brought 
about  by  the  discrete  elements.  Other  conclusions  of  the  study  are:  (i)  the  HHC  is  effective 
in  quickly  reducing  vibrations,  at  least  at  its  design  condition;  (ii)  a  linearized  model  of 
helicopter  dynamics  is  adequate  for  HHC  design,  as  long  as  the  periodicity  of  the  system  is 
correctly  taken  into  account,  i.e.,  periodicity  is  more  important  than  nonlinearity,  at  least  for 
the  mathematical  model  used  in  this  study;  and  (iii)  when  discrete  and  continuous  systems 
are  both  stable,  the  predicted  HHC  control  harmonics  are  in  good  agreement,  although  the 
initial  transient  behavior  can  be  considerably  different. 
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Chapter  2 


Key  accomplishments 

2.1  Research  accomplishments 

The  key  accomplishments  of  the  project  are  briefly  summarized  below: 

1.  First  systematic  study  of  the  effect  of  zeros  in  rotorcraft  aeromechanics.  The  “zeros” 
of  a  dynamic  system  play  a  key  role  in  the  closed-loop  behavior  of  the  system.  For  ex¬ 
ample,  if  one  or  more  of  the  zeros  have  positive  real  parts,  a  control  system  designed  to 
stabilize  the  system  could  instead  make  it  unstable.  The  rotorcraft  dynamics  commu¬ 
nity  had  been  able  to  mostly  neglect  the  issue  because  rotor  systems  were  not  actively 
controlled.  As  rotor  active  controls  become  technologically  feasible  and  desirable,  the 
calculation  and  study  of  zeros  must  become  a  routine  step  of  rotor  design. 

2.  First  study  of  the  effects  of  closed-loop  HHC/IBC  on  the  aeroelastic  stability  of  a  heli¬ 
copter  rotor.  Despite  over  three  decades  of  research  in  HHC/IBC,  no  information  was 
available  on  how  an  HHC/IBC  system  would  affect,  for  example,  the  in-plane  damping 
of  a  helicopter  rotor.  No  results  showing  the  effects  of  HHC/IBC  on  stability  eigenval¬ 
ues  or  Floquet  characteristic  exponents  for  rotor  and/or  fuselage  modes  had  ever  been 
presented.  This  is  clearly  a  vital  piece  of  information  for  the  future  use  of  active  rotor 
controls  in  configurations  with  inherently  low  inplane  damping  such  as  hingeless  and 
bearingless  rotors. 
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3.  Presentation  to  the  rotorcraft  community  of  techniques  for  the  analysis  of  systems 
with  periodic  coefficients  developed  in  other  areas  of  engineering.  Several  time-  and 
frequency-domain  techniques,  developed  in  the  last  two  decades  primarily  in  the  signal 
processing  and  control  engineering  communities,  have  been  introduced  to  the  rotorcraft 
community  through  relevant,  helicopter-oriented  application  examples. 

4.  First  study  of  the  effects  of  closed-loop  HHC/IBC  on  the  aeroelastic  stability  of  a  he¬ 
licopter,  considered  as  a  discrete  system,  A  typical  model  of  a  helicopter  with  HHC 
or  IBC  contains  continuous  portions  (e.g.,  rotor  and  fuselage  dynamics)  and  discrete 
portions  (e.g.,  harmonic  analysis  to  extract  vibratory  components,  HHC  inputs  up¬ 
dated  once  per  revolution,  etc.).  The  stability  study  mentioned  as  item  3.  above  was 
carried  out  by  approximately  modeling  some  discrete  elements  as  continuous  and  sim¬ 
ply  neglecting  the  rest.  Subsequently,  we  have  modeled  the  entire  system  as  discrete. 
As  a  consequence,  we  have  developed  the  first  rigorous  solution  of  the  aeroelastic  sta¬ 
bility  problem  for  a  helicopter  with  HHC  or  IBC,  including  the  effects  of  the  discrete 
elements,  and  of  different  sampling  and  update  rates  in  the  system.  Also,  this  was 
essentially  the  first  solution  of  the  aeroelastic  equations  as  discrete  systems  (a  previous 
example  in  the  literature  was  much  more  limited  in  scope). 

In  all  the  accomplishments  just  listed,  the  mathematical  model  of  the  helicopter  was 
of  realistic  complexity,  and  not  an  idealized,  simplified  representation.  In  most  cases,  the 
model  contained  more  than  35  degrees  of  freedom,  including  finite  element-based  rotor  blade 
coupled  flag- lag-torsional  dynamics,  full  six  degree-of-freedom  rigid  body  dynamics,  and 
inflow  dynamics. 

Finally,  Dr.  Rendy  Cheng,  as  part  of  his  Ph.D.  dissertation  at  the  University  of  Maryland 
(with  Dr.  Ccli  as  academic  advisor),  and  in  collaboration  with  Dr.  Mark  Tischler  of  the  U.  S. 
Army  Aeroflightdynamics  Directorate,  Ames  Research  Center,  developed  a  model  for  the 
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closed-loop  analysis  of  a  UH-60  with  an  IBC  system.  This  work,  intended  to  support  the 
RADICL  project  at  Ames,  was  not  directly  part  of  the  ARO  sponsored  research.  However, 
the  results  of  the  ARO  research  helped  clarify  the  results  of  Dr.  Cheng’s  research  from  a 
theoretical  point  of  view. 

2.2  Publications 

The  research  performed  in  this  project  resulted  in  the  following  conference  and  journal 
publications: 

2.2.1  Conference  publications 

1.  Lovera,  M.,  Colaneri,  P.,  Celi,  R.,  and  Bittanti,  S.,  “On  the  Role  of  Zeros  in  Rotorcraft 
Aeromechanics,”  Proceedings  of  the  58th  Annual  Forum  of  the  American  Helicopter 
Society ,  Montreal,  Canada,  June  2002. 

2.  Lovera,  M.,  Colaneri,  P.,  and  Celi,  R.,  “Periodic  Control  Issues  in  the  Higher  Harmonic 
Control  of  Helicopter  Rotors,”  Proceedings  of  the  2003  American  Control  Conference, 
Denver,  USA,  June  2003. 

3.  Lovera,  M.,  Colaneri,  P.,  Malpica,  C.,  and  Celi,  R.,  “Closed-Loop  Aeromechanical 
Stability  Analysis  for  a  Hingeless  Rotor  Helicopter  With  HHC  or  IBC,”  Proceedings  of 
the  29th  European  Rotorcraft  Forum,  Friedrichshafen,  Germany,  September  2003. 

4.  Colaneri,  P.,  Celi,  R.,  and  Bittanti,  S.,  “Constant-Coefficient  Representations  of  Periodic- 
Coefficient  Discrete  Linear  Systems,”  Proceedings  of  the  AHS  fth  Decennium  Specialist 
Conference  on  Aeromechanics,  San  Francisco,  CA,  January  2004. 

5.  Lovera,  M.,  Colaneri,  P.,  Malpica,  C.,  and  Celi,  R.,  “Discrete-Time,  Closed-Loop 
Aeromechanical  Stability  Analysis  of  Helicopters  With  Higher  Harmonic  Control,”  Pro- 
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ceedings  of  the  60th  Annual  Forum  of  the  American  Helicopter  Society,  Baltimore,  MD, 
June  2004. 

2.2.2  Journal  publications 

1.  Lovera,  M.,  Colaneri,  P.,  and  Celi,  R.,  “On  the  Role  of  Zeros  in  Rotorcraft  Aerome¬ 
chanics,”  to  appear  in  the  October  2004  issue  of  the  Journal  of  the  American  Helicopter 
Society. 

2.  Lovera,  M.,  Colaneri,  P.,  Malpica,  C.,  and  Celi,  R.,  “Closed-Loop  Aeromechanical  Sta¬ 
bility  of  a  Hingeless  Rotor  Helicopter  with  Higher  Harmonic  Control,”  accepted  for 
publication  in  the  Journal  of  Guidance,  Control,  and  Dynamics. 

3.  Colaneri,  P.,  Celi,  R.,  andBittanti,  S.,  “Constant-Coefficient  Representations  of  Periodic- 
Coefficient  Discrete  Linear  Systems,”  to  be  submitted  for  publication  in  the  Journal 
of  the  American  Helicopter  Society. 

4.  Lovera,  M.,  Colaneri,  P.,  Malpica,  C.,  and  Celi,  R.,  “Discrete-Time,  Closed-Loop 
Aeromechanical  Stability  Analysis  of  Helicopters  With  Higher  Harmonic  Control,”  to 
be  submitted  for  publication  in  the  Journal  of  Guidance,  Control,  and  Dynamics. 
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Chapter  3 


On  the  Role  of  Zeros  in  Rotorcraft  Aeromechanics 

3.1  Introduction  and  Problem  Statement 

Various  types  of  active  control  of  a  helicopter  rotor  have  been  proposed  as  a  means  to 
improve  main  rotor  performance,  and  reduce  vibratory  loads  and  noise  (see  e.g.,[20,  37]). 
The  rotor  controls  would  be  based  on  nonrotating  actuators  mounted  in  the  fuselage  or  on 
rotating  actuators,  either  mounted  on  the  rotating  portion  of  the  swashplate  and  replacing 
the  conventional  rigid  pitch  links,  or  on  the  rotor  blades  themselves  in  the  form  of  trailing 
edge  flaps,  spoilers,  elevons,  or  active  tips. 

Traditionally,  the  dynamics  community  has  been  primarily  interested  in  the  open  loop 
behavior  of  the  rotor.  For  example,  the  behavior  of  the  stability  eigenvalues  (or  of  the  Flo- 
quet  characteristic  exponents,  if  the  periodicity  of  the  rotor  equations  of  motion  is  taken 
into  account)  has  been  studied  extensively,  and  it  is  now  well  known  how  they  are  influ¬ 
enced  by  flight  condition  and  rotor  configuration  [19,  10].  Stability  eigenvalues  and  Floquet 
characteristic  exponents  can  be  seen  as  the  poles  of  the  system. 

When  a  feedback  loop  is  closed  around  the  rotor,  not  only  the  poles,  but  also  the  zeros 
of  the  system  become  significant.  For  example,  zeros  in  the  right-half-plane  (RHP),  or 
nonminimum  phase  (NMP)  zeros,  can  limit  the  performance  achievable  through  feedback 
control,  and  more  or  less  precise  pole-zero  cancellation  can  increase  or  decrease  intermodal 
coupling.  For  a  single-input  single-output  (SISO)  system,  high  values  of  the  feedback  gains 
are  desirable  for  good  tracking  of  the  desired  response,  and  to  reduce  the  sensitivity  to 
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actuator  noise,  modeling  uncertainties,  and  changes  in  operating  conditions.  As  the  gains 
increase,  the  zeros  of  the  system  tend  to  “attract”  some  of  the  closed-loop  poles  of  the 
systems,  and  if  the  system  has  unstable,  or  NMP,  zeros,  closed-loop  instabilities  can  result. 
Many  of  these  concepts  carry  over  to  the  multi-input,  multi-output  (MIMO)  case. 

Clearly,  a  fundamental  understanding  of  coupled  rotor-fuselage  zeros  is  a  key  ingredient 
for  the  success  of  closed-loop,  active  rotor  control.  Yet,  with  the  partial  exception  of  Refs.  [24] 
and  [5],  no  systematic  study  on  this  topic  has  appeared  in  the  literature.  Ref.  [5]  used  a 
relatively  simple,  isolated  rotor  model,  and  showed  that  NMP  zeros  can  indeed  be  present, 
both  in  hover  and  in  forward  flight. 

A  rigorous  study  of  the  aeroclastic  stability  of  a  rotor  system  in  forward  flight  should 
include  the  effects  of  periodic  coefficients  in  the  governing  equations.  A  simpler,  constant 
coefficient  approximation  can  be  obtained  by  performing  a  Multiblade  Coordinate  Trans¬ 
formation  (MCT)  and  neglecting  the  residual  periodicity  of  the  equations:  the  poles,  or 
stability  eigenvalues,  of  this  approximate  model  are  typically  accurate  up  to  an  advance  ra¬ 
tio  of  /i  =  0.3  —  0.35.  No  studies  have  been  published  on  whether  such  approximations  are 
also  accurate  for  the  zeros ,  and,  if  this  is  the  case,  up  to  what  advance  ratio. 

In  the  literature  related  to  the  analysis  of  systems  with  periodic  coefficients  the  problem 
of  defining  and  computing  system  zeros  has  been  studied  only  in  recent  years.  The  first 
contributions  were  devoted  to  the  case  of  discrete  time  systems  (see  [7],  [23]  and  [15]).  On 
the  other  hand,  the  case  of  continuous  time  systems  has  been  investigated  for  the  first  time 
in  Ref.  [16].  In  this  paper,  a  system  theoretic  definition  for  the  zeros  of  a  continuous  time 
linear  system  with  periodic  coefficients  is  given. 

In  light  of  the  preceding  discussion,  the  objectives  of  the  present  chapter  are: 

1.  To  summarize  the  main  properties  of  the  zeros  of  linear  systems  with  periodic  coeffi¬ 
cients,  especially  as  they  apply  to  a  helicopter  rotor  system,  and  their  implication  on 
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rotor  design; 


2.  To  present  methods  for  the  calculation  of  the  zeros,  both  for  governing  equations  with 
constant  and  with  periodic  coefficients,  and  to  discuss  the  accuracy  of  approximate, 
constant  coefficient  models;  and 

3.  To  present  results  for  realistic  hingeless  rotor  configurations,  both  in  hover  and  forward 
flight. 

A  fundamental  difference  between  the  calculation  of  poles  and  of  zeros  is  that  the  former 
is  only  based  on  the  homogeneous  system,  whereas  specific  outputs  and  controls  need  to  be 
defined  for  the  latter.  Therefore,  any  matrix  of  test  cases  is  multidimensional  and  potentially 
very  large.  The  results  of  this  chapter  will  be  limited  to  rotating  swashplate  inputs,  and  the 
output  will  be  limited  to  fuselage  and  rotor  states.  Other  possible  outputs,  such  as  hub  loads 
or  fuselage  vibrations  at  points  other  than  the  center  of  gravity,  will  not  be  considered. 

3.2  Helicopter  simulation  model 

The  baseline  simulation  model  used  in  this  chapter  is  a  nonreal-time,  blade  element  type, 
coupled  rotor-fuselage  simulation  model  [38].  The  fuselage  is  assumed  to  be  rigid  and  dy¬ 
namically  coupled  with  the  rotor.  A  total  of  nine  states  describe  fuselage  motion  through 
nonlinear  Euler  equations.  Fuselage  and  blade  aerodynamics  are  described  through  tables 
of  aerodynamic  coefficients,  and  no  small  angle  assumption  is  required.  A  coupled  flap-lag- 
torsion  clastic  rotor  model  is  used.  Blades  are  modeled  as  Bernoulli-Euler  beams.  The  rotor 
is  discretized  using  finite  elements,  with  a  modal  coordinate  transformation  to  reduce  the 
number  of  degrees  of  freedom.  The  elastic  deflections  are  not  required  to  be  small.  Blade 
element  theory  is  used  to  obtain  the  aerodynamic  characteristics  on  each  blade  section. 
Quasi-steady  aerodynamics  is  used,  with  a  3-state  dynamic  inflow  model. 


21 


The  trim  procedure  is  the  same  as  in  Ref.  [8].  The  rotor  equations  of  motion  are  trans¬ 
formed  into  a  system  of  nonlinear  algebraic  equations  using  a  Galerkin  method.  The  algebraic 
equations  enforcing  force  and  moment  equilibrium,  the  Euler  kinematic  equations,  the  inflow 
equations  and  the  rotor  equations  are  combined  in  a  single  coupled  system.  The  solution 
yields  the  harmonics  of  a  Fourier  expansion  of  the  rotor  degrees  of  freedom,  the  pitch  control 
settings,  trim  attitudes  and  rates  of  the  entire  helicopter,  and  main  and  tail  rotor  inflow. 
Linearized  models  are  extracted  numerically,  by  perturbing  rotor,  fuselage,  and  inflow  states 
about  a  trimmed  equilibrium  position. 

3.3  Zeros  of  linear  time-periodic  systems 

Consider  the  continuous-time  multivariable  linear  periodic  system 

x(t)  =  A(t)x(t)  +  B(t)u(t)  ,  x 

y{t)  —  C(t)x(t)  +  D(t)u(t) 

where  all  the  matrices  have  a  common  fundamental  period  T,  such  that  it  is,  for  example 
A(t)  =  A(t  +  T),  and  the  state  vector  x(t),  the  control  vector  u(t),  and  the  output  vector 
y(t),  have  dimension  n,  m,  and  p,  respectively. 

Introduce  then  the  class  of  Exponentially  Modulated  Periodic  (EMP)  signals  [41].  The 
(complex)  signal  u(t)  is  said  to  be  EMP  of  period  T  and  modulation  s  if 

u(t)  =  UneSnt  = eSt  UnejnQt  (3-2) 

n=0, ±1, ±2,...  n=0,±l,±2,... 

where  t  >  0,  sn  —  s  +  jnfl,  and  s  is  a  complex  scalar. 

The  class  of  EMP  signals  is  a  generalization  of  the  class  of  T-periodic  signals  (i.e.,  signals 
with  period  T)  in  the  sense  that  u(t)  can  be  written  as  u(t)  =  u[t)est  where  u(t)  is  a  T- 
periodic  signal.  Consequently,  an  EMP  signal  with  s  =  0  is  just  an  ordinary  time-periodic 
signal.  In  much  the  same  way  as  a  time  invariant  system  subject  to  a  (complex)  exponential 
input  has  an  exponential  steady-state  response,  a  periodic  system  subject  to  an  EMP  input 
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has  an  EMP  steady-state  response.  In  such  a  response,  all  signals  of  interest  (x,  x,  y )  can 
be  expanded  as  EMP  signals.  By  deriving  Fourier  expansions  for  A(t),  B(t),  C(t)  and  D(t), 
it  is  possible  to  prove  that  the  EMP  steady-state  response  of  the  system  can  be  expressed 
as  an  infinite  dimensional  matrix  equation  with  constant  elements  (see  Ref.  [41]  for  details). 

This  chapter  will  focus  on  periodic  systems  with  equal  number  of  inputs  and  outputs,  i.e., 
square  systems  (m  =  p).  Moreover,  the  assumption  will  be  made  that  the  given  system  is 
reachable  and  observable  (a  formal  definition  of  the  structural  properties  of  periodic  systems 
can  be  found  in  Ref.  [2]).  With  this  last  assumption  it  is  not  necessary  to  distinguish  between 
transmission  and  invariant  zeros,  and  therefore  it  is  possible  to  refer  simply  to  the“  zeros  of 
a  periodic  system”  without  any  additional  specification.  A  complex  number  z  is  said  to  be  a 
zero  of  the  system  of  Eq.  (3.1)  if  there  exist  two  T-periodic  vectors  x(t)  and  u(t)  (not  both 
identically  zero)  such  that 

'(a  +  z)I-A(t)  ~B(t)  1  [  x(t)  1  / oox 

C{t )  D{t)  \  [  u(t )  J  1  j 

where  a  is  the  derivative  operator  d/dt.  Equation  (3.3)  can  be  rewritten  in  terms  of  the 
so-called  blocking  property  condition.  Precisely,  the  complex  number  z  is  a  zero  of  the 
system  of  Eq.  (3.1)  if  and  only  if  there  exist  an  exponentially  modulated  periodic  signal 
u(t)  =  ^n=0±i±2  uneSnt  and  an  initial  state  x(0)  such  that  the  corresponding  output  y(t) 
is  identically  zero. 

The  above  definition  of  periodic  zeros,  based  on  the  blocking  property,  is  not  the  most 
appropriate  one  when  dealing  with  the  problem  of  computing  zeros.  Therefore,  a  different 
definition  that  leads  to  a  viable  computational  approach  will  be  discussed  next. 

First,  introduce  the  operator 

LAh{t)  =  h(t)A(t)  +  ^rh(t)  (3.4) 

dt 

where  A(t)  is  a  square  n  x  n  matrix,  and  h(t)  is  a  vector  function  of  compatible  dimensions. 
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Then,  with  reference  to  the  system  of  Eq.  (3.1)  consider  the  following  definition  of  “relative 
degree  r(t)n  [16]. 

For  i  =  1,...  ,p,  suppose  that  there  exists  an  integer  valued  function  rt  =  rjt),  1  < 
ri(t)  <  n  such  that,  indicating  by  cjt)  the  i-th  row  of  the  matrix  C(t): 

1.  Ci{r)B{r )  =  L,ACi(r)B(T )  =  . . .  =  Lrj^2Ci{T)B{r)  =  0  for  all  r  in  a  neighborhood  of  t. 

2.  in  every  neighborhood  of  t  there  exists  r  such  that  L([_1Cj(r)i?(r)  ^  0. 


Then,  the  system  is  said  to  have  relative  degree  r(t )  =  [ri(t)  . . .  rp(t)]  at  time  t. 

The  above  definition  of  relative  degree  is  the  natural  extension  to  time-periodic  systems 
of  the  same  concept  for  time-invariant  systems.  Indeed,  for  a  SISO  LTI  system  of  order  n 
described  by  a  transfer  function 

n(  ,  N(s)  n0sm  +  nism~1  +  . . .  +  nm 

G(s)  =  m =  s» +</„•.-! +...+d„  (3-5) 

with  n0  ^  0,  the  relative  degree  rLTi,  defined  as 

r  lti  —  n  —  m  (3.6) 


(i.e.,  equal  to  the  difference  between  the  degree  of  the  polynomial  at  the  denominator  of  the 
transfer  function  and  that  at  the  numerator  or,  equivalently,  to  the  difference  between  the 
number  of  poles  and  the  number  of  zeros)  can  be  given  a  state  space  interpretation  which 
coincides  with  the  given  definition.  The  transfer  function  G(s),  Eq.  (3.5),  can  be  written  as 


G(s) 


N(s) 
D(s ) 


C{sl  -  A)~1B  +  D 


N(s) 


D(s) 


(3.7) 


where  N(s)/D(s)  '=  C(sl  —  A )  1B,  and  N(s )  and  D(s)  are,  respectively,  polynomials  in 
s  of  degree  n  and  n  —  1.  If  the  D  matrix  (in  the  SISO  case,  a  1  by  1  matrix)  is  not  equal 
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to  zero,  then  N(s)  =  N(s)  +  D(s)D  is  a  polynomial  of  degree  n,  and  the  relative  degree  is 
r  —  0.  If  D  =  0  but  CB  0  it  can  be  shown  (see,  for  example,  Ref.  [34])  that  the  relative 
degree  is  r  =  1.  If  D  —  0  and  CB  =  0,  but  CAB  ^  0,  the  relative  degree  is  r  =  2,  and  so 
on.  For  a  MIMO  system  the  relative  degree  becomes  a  vector,  each  element  of  which  gives 
the  relative  degree  of  each  output  channel. 

In  principle,  for  a  generic  periodic  system  the  relative  degree  is  a  function  of  time.  This 
means  that  in  general  the  number  of  zeros  of  the  system  might  be  itself  a  function  of  time. 
However,  in  most  practical  problems  the  number  of  zeros  is  fixed,  the  relative  degree  r 
is  independent  from  t,  and  r  can  be  calculated  using  the  same  tests  on  A,  B ,  C,  and  D 
previously  mentioned;  the  only  condition  is  that  if  any  of  A,  B,  C ,  D  becomes  equal  to  zero, 
it  does  so  only  for  isolated  instants  in  time,  and  not  time  intervals  (see  Ref.  [16]  for  a  more 
mathematically  rigorous  discussion).  This  leads  to  the  notion  of  uniform,  relative  degree. 

Assuming  first  that  the  system  of  Eq.  (3.1)  has  relative  degree  r(t),  one  can  define  an 
instantaneous  interaction  m,atrix  T(t)  of  the  system  as 

r(t)  =  :  (3,8) 

Then  the  system  of  Eq.  (3.1)  is  said  to  have  uniform  relative  degree  r  =  [rq  . . .  rp]  if  r(t)  =  r, 
1  <  Ti  <  n,  i  =  1, ...  ,p,  the  rank  of  the  matrix  T(t)  is  equal  to  p  for  every  value  of  t,  and 
r (t)  has  a  bounded  inverse.  It  can  be  shown  (see  Ref.  [16])  that,  for  such  a  system,  the  zeros 
coincide  with  n  —  r  periodic  poles  of  the  system 

x(t)  =  R(t)x(t )  (3.9) 

where  the  matrix  R(t)  is  defined  as 

R(t)  =  Aft)  -  B^Y^Qit)  (3.10) 
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and  the  p  by  n  matrix  Q(t)  is  given  by 

'  LTcAt)  " 

Qit)  =  :  ,  (3.11) 

.  LrX~lcP(t)  _ 

the  remaining  r  poles  being  zero.  As  in  the  case  of  the  definition  of  periodic  poles,  pe¬ 
riodic  zeros  can  be  described  both  in  terms  of  characteristic  multipliers  and  characteristic 
exponents. 

The  previously  introduced  definitions  related  to  periodic  zeros  hold  for  a  generic  square 
continuous  time  linear  periodic  system,  i.e.,  a  system  with  equal  number  of  inputs  and 
outputs.  The  extension  of  the  above  definition  to  nonsquare  periodic  systems  is  presently  a 
subject  of  research. 

3.4  Closed  loop  zeros  of  periodic  systems 

Another  interesting  issue  in  the  study  of  periodic  zeros  is  their  role  in  determining  the  closed 
loop  behavior  of  control  systems  in  high  gain  situations,  ft  turns  out  that  periodic  zeros  have 
the  same  property  of  their  time-invariant  counterpart  of  attracting  closed  loop  poles  for  high 
gain.  In  fact,  it  can  be  shown  that  a  finite  zero  of  a  square  periodic  system  can  be  seen  as 
the  limit  of  a  closed  loop  pole  as  the  gain  matrix  goes  to  infinity. 

This  result  plays  a  fundamental  role  in  explaining  why  the  closed  loop  behavior  of  an 
”  almost  time  invariant”  LTP  system  can  turn  out  to  be  completely  different  from  the  one 
of  its  LTI  approximation.  To  illustrate  this  fact,  consider  the  following  simple  example:  a 
linear  time  periodic  system  is  given,  which  is  described  by  the  matrices 

A(t)=\  -1+si,;(‘>  0  lS(t)=f-1-cosWl  (3.12) 

1  —  cos  (t)  — 3  J  L  2_  smp) 

C{t)  =[01]  ,D(t)  =  0.  (3.13) 
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As  A{t)  is  lower  triangular,  it  is  easy  to  see  that  the  poles  (characteristic  exponents)  of  the 
periodic  system  are  given  by  —1  and  —3  and  coincide  with  the  poles  of  the  corresponding 
time-invariant  approximation  of  the  system. 

Consider  now  the  closed  loop  system 


A(t)  =  A{t)  +  kB(t)C(t) 

—  1  +  sin(f)  — k  —  kcos(t) 

1  —  cos (t)  —3  +  2 k  —  k  sin (t) 

and  its  time-invariant  approximation 


Aq  — 


-1  —k 
1  —3  +  2k 


(3.14) 


(3.15) 


The  constant  coefficient  A0  is  asymptotically  stable  for  k  <  2,  however  the  time-periodic 
A(t)  is  asymptotically  stable  only  for  k  <  1.6,  so  in  this  case  the  closed  loop  stability 
analysis  based  on  the  LT1  approximation  can  be  misleading,  even  if  the  open  loop  periodic 
poles  coincide  with  the  time-invariant  ones.  As  the  two  systems  have  the  same  poles,  the 
discrepancy  in  the  closed  loop  behavior  can  only  be  interpreted  by  looking  at  the  zeros.  The 
LT1  system  has  relative  degree  equal  to  1  (transfer  function  with  two  poles  and  one  zero 
or,  equivalently,  D0  =  0  and  Co-Bo  7^  0),  and  the  same  holds  for  the  time-periodic  one,  as 
the  product  C(t)B(t )  is  generically  nonzero,  as  required  by  the  definition  of  relative  degree 
r(t)  previously  introduced.  The  zero  of  the  time  invariant  system  [A0,  B0 ,  Co,  D0\  is  given  by 
zlti  =  —0.5,  while  the  zero  of  the  time  periodic  system  [A(t),B(t),C(t),D(t)],  computed 
according  to  the  definitions  in  the  previous  Section,  is  given  by  zltp  =  —0.69. 

Clearly,  the  different  behavior  of  the  closed  loop  system  in  the  two  cases  is  clue  to  the 
different  structure  of  the  zero  dynamics,  so  in  this  example  a  correct  analysis  of  the  closed 
loop  behavior  of  the  system  can  be  carried  out  only  by  taking  into  account  the  periodicity 
of  the  dynamics.  To  this  purpose,  the  following  Section  will  show  how  to  compute  periodic 
zeros  for  the  case  of  rotor-fuselage  models. 
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3.5  Calculation  of  coupled  rotor-fuselage  zeros 


The  linearized  coupled  rotor-fuselage  system  can  be  modeled  as  a  continuous-time  linear 
periodic  system  of  the  type  of  Eq.  (3.1),  where  the  matrices  [A(t),  B(t),C(t),  D(t)\  have  a 
common  fundamental  period  T,  equal  to  1/N  rotor  revolutions.  While,  for  a  given  flight 
condition  and  helicopter  configuration,  the  matrix  A(t)  is  fixed,  matrices  B(t),C(t),D(t) 
vary  according  to  the  input-output  configuration  one  wants  to  study.  For  simplicity,  this 
chapter  will  concentrate  on  the  analysis  of  SISO  rotor  fuselage  models,  however  the  extension 
of  the  results  presented  here  to  the  more  general  case  of  a  (square)  MIMO  model  does  not 
pose  any  theoretical  or  computational  problems. 

As  described  in  the  previous  Section,  the  simplest  way  to  compute  the  zeros  of  a  periodic 
system  is  to  cast  the  problem  into  that  of  computing  of  the  characteristic  exponents  of  an 
associated  periodic  matrix  R(t).  While  the  general  form  of  R(t)  given  in  equation  (3.10)  is 
a  complicated  one,  the  problem  is  considerably  simpler  if  one  considers  the  particular  cases 
of  systems  with  uniform  relative  degree  r  —  0  ( D(t )  generically  non  singular)  and  r  —  1 
(D(t)  =  0  for  all  t,  and  C(t)B(t )  generically  non  singular). 

If  D(t)  is  invertible  for  each  t,  then  the  zeros  of  the  system  coincide  with  the  characteristic 
exponents  of  the  inverse  system,  i.e.,  of  matrix 

R(t)  =  A(t)  -  B(t)D(t)~lC(t)  (3.16) 


If  the  system  has  a  singular  D(t)  ( D(t )  =  0  in  the  SISO  case)  and  C(t)B(t )  is  generically 
non  singular  ( C(t)B(t )  ^  0  for  all  t  in  the  SISO  case),  then  the  zeros  of  the  system  coincide 
with  the  characteristic  exponents  of 


R(t)  =  A(t)  -  B(t)[C(t)B(t)]^ 


dCjt) 

dt 


+  C(t)A(t) 


(3.17) 


except  for  one  at  the  origin.  For  the  determination  of  r,  it  has  been  assumed  that  the  relative 
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degree  of  the  periodic  models  coincides  with  the  one  of  the  corresponding  LTI  approxima¬ 
tions,  which  can  be  computed  very  easily. 

The  problem  of  computing  the  characteristic  multipliers  of  a  linear  periodic  system  has 
been  extensively  studied  in  the  literature,  and  a  number  of  approaches  have  been  proposed. 
In  particular,  four  classes  of  computational  methods  can  be  defined  (see  Ref.  [30]  for  a 
detailed  overview): 

•  Direct  integration  methods  (single  or  multiple  shooting  [24]). 

•  Frequency  domain  methods  [41,  43]. 

•  Symbolic  methods  [36]. 


•  Linear  algebra  methods  (periodic  Schur  decomposition  [6,  31]). 

From  a  practical  viewpoint,  as  the  helicopter  model  used  in  this  chapter  has  53  state 
variables,  the  use  of  the  symbolic  and  frequency  response  methods  is  probably  not  a  viable 
option.  Therefore,  only  the  direct  integration  approach  and  the  method  based  on  the  periodic 
Schur  decomposition  have  been  used.  The  main  difference  between  the  two  approaches  can 
be  best  understood  by  looking  at  the  operations  needed  to  determine  the  characteristic 
multipliers.  The  multipliers  are  defined  as  the  eigenvalues  of  the  monodromy  matrix  (i.e.,  the 
Floquet  transition  matrix  at  the  end  of  one  period)  of  the  periodic  system  under  investigation. 
Therefore,  the  conventional  (direct  integration)  approach  to  the  problem  consists  of  the 
following  steps: 


1.  Solving  the  initial  value  problem 

0)  =  0) 

3*r(0,  0)  =  In 

over  the  interval  from  0  to  T,  to  compute  the  monodromy  matrix  M 
numerical  integration  algorithm  can  be  used  to  this  purpose. 


(3.18) 
<hR(T,0).  Any 
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2.  Computing  the  characteristic  multipliers  as  the  eigenvalues  of  M . 


This  scheme  has  been  successfully  used  for  the  computation  of  poles  of  rotorcraft  dy¬ 
namics,  as  the  open  loop  dynamics  do  not  give  rise  to  any  major  numerical  problems  in  the 
solution  of  Eq.  (3.18)  (the  coupled  rotor-fuselage  dynamics  is  at  most  only  mildly  unstable 
and  not  particularly  stiff).  Unfortunately,  the  situation  can  be  considerably  different  when 
dealing  with  the  problem  of  computing  zeros.  In  fact,  zeros  that  are  very  large  and  pos¬ 
itive  can  often  occur.  In  this  case  the  zero  dynamics  will  be  strongly  nonminimum  phase 
(i.e.,  unstable)  and  stiff,  leading  to  potentially  severe  numerical  problems  in  the  solution  of 
Eq.  (3.18).  For  example,  if  the  real  part  of  one  or  more  nonminimum  phase  zeros  is  suffi¬ 
ciently  large,  the  integration  required  to  compute  the  monodromy  matrix  will  fail  because 
of  numerical  overflow  before  reaching  the  end  of  the  period.  Also,  the  monodromy  matrix 
associated  with  the  zero  dynamics  can  be  poorly  scaled.  These  problems  do  not  necessarily 
occur  in  every  single  zero  calculation,  but  they  can  easily  occur  in  otherwise  well  behaved, 
practical  problems.  As  a  consequence,  it  is  desirable  to  use  a  computational  scheme  that  can 
deal  with  these  situations,  should  they  occur.  The  so-called  “periodic”  Schur  decomposition 
is  one  such  scheme,  based  on  the  idea  of  avoiding  the  explicit  calculation  of  the  monodromy 
matrix  M.  The  main  idea  behind  this  approach  is  the  following. 

Consider  a  sequence  of  time  instants  4,  k  —  1, . . .  ,  K ,  with  t0  =  0,  4+i  >  4,  and  tx  =  T, 
and  let 

Mk=f  $R(tk,tk.  i)  k  =  l,...,K  (3.19) 

Then  clearly  one  has 

M  =  M KMK~\  . . .  M\  (3.20) 

If  each  of  the  Mk  matrices  is  diagonal  or  lower  (upper)  triangular,  then  the  characteristic 
multipliers  are  given  by  the  product  of  the  eigenvalues  of  the  Mk  matrices.  Because  in 
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general  this  is  not  the  case,  the  periodic  Schur  decomposition  algorithm  can  be  used  to 
apply  a  suitable  sequence  of  unitary  operations  to  the  Mk  matrices  in  order  to  transform 
each  of  them  to  upper  triangular  form  (except  one  which  can  be  in  quasi  upper  triangular 
form).  It  can  be  proved  (see  Refs. [6,  31,  21]  for  details)  that  this  decomposition  of  the  Mk 
matrices  corresponds  to  an  implicit  Schur  decomposition  for  M,  from  which  the  characteristic 
multipliers  can  be  easily  be  computed.  Therefore,  the  overall  numerical  scheme  for  the 
computation  of  the  characteristic  multipliers  proceeds  as  follows: 

1.  Select  K  and  the  time  instants  tk  (usually  equally  spaced)  and  compute  the  K  tran¬ 
sition  matrices  from  tk~ i  to  tk.  Note  that  if  K  is  chosen  sufficiently  large,  then  the 
assumption  that  R(t)  is  approximately  constant  over  each  interval  can  be  introduced 
and  each  of  the  Mk  can  be  computed  as 

Mk  «  e*(tk-i)[*k-tk-i]  (3.21) 

2.  Use  the  periodic  Schur  decomposition  algorithm  to  determine  the  upper  triangular  (or 
quasi  upper  triangular)  form  of  each  of  the  Mk  matrices  in  the  appropriate  order. 

3.  Compute  the  eigenvalues  of  the  monodromy  matrix  M  from  the  triangular  form  of  the 
matrices  Mk. 

This  computational  scheme  has  two  major  advantages  over  the  direct  integration  method: 
firstly,  it  does  not  require  the  numerical  integration  of  the  dynamics,  which  can  be  a  critical 
step  when  dealing  with  zero  dynamics  (stiffness);  secondly,  it  avoids  the  explicit  formation  of 
M  and  applies  (almost)  only  unitary  operations  to  the  data,  thus  minimizing  the  propagation 
of  numerical  errors.  The  results  presented  in  the  following  Section  have  been  obtained  using 
this  computational  method,  in  particular  by  resorting  to  the  implementation  of  the  periodic 
Schur  algorithm  described  in  [31]. 
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3.6  Results 


All  the  results  presented  in  this  section  refer  to  a  hingeless  rotor  configuration  similar  to  the 
Eurocopter  BO-105.  The  linearized  system  has  53  states:  9  for  the  rigid  body  motion;  40 
to  model  five  elastic  coupled  modes,  for  each  of  the  four  rotor  blades;  3  for  the  main  rotor 
dynamic  inflow  model;  and  one  for  the  tail  rotor  dynamic  inflow  model.  In  this  chapter,  an 
analysis  of  the  complex  plane  representations  of  the  poles  and  zeros  of  the  SISO  transfer 
function  from  collective  input  to  the  acceleration  w  of  the  center  of  gravity  of  the  helicopter 
is  presented.  This  is  one  of  the  open  loop  transfer  functions  that  would  have  to  be  taken 
into  account  in  the  implementation  of  IBC  for  vibration  reduction. 

Three  conhgnrations  will  be  analyzed,  namely: 

1.  A  baseline  configuration,  with  value  of  the  thrust  coefficient  Ct/ct  =  0.07  and  the 
basic  soft-in-plane  rotor  with  a  fundamental  rotating  lag  frequency  of  approximately 
0.7  /rev; 

2.  A  heavier  configuration,  with  Ct/ct  =  0.10,  and  the  same  soft-in-plane  rotor;  and 

3.  A  stiff- in-plane  configuration,  with  the  baseline  value  of  Ct/ct  =  0.07  and  a  fundamen¬ 
tal  rotating  lag  frequency  of  approximately  1.4/rev.  This  configuration  has  often  been 
used  in  rotary-wing  aeroelasticity  studies  because  it  becomes  unstable  at  high  advance 
ratios  [18]. 

The  results  will  cover  a  speed  range  from  hover  to  up  to  150  kts,  depending  on  the  configu¬ 
ration.  The  upper  speed  limit  is  determined  in  each  case  by  the  maximum  speed  at  which  a 
converged  trim  solution  can  be  achieved.  Because  a  value  of  /j  =  0.1  corresponds  to  a  speed 
of  about  42  kts,  the  maximum  advance  ratio  considered  in  the  results  is  just  above  /c  =  0.35. 

All  the  states  are  defined  in  a  nonrotating  coordinate  system.  The  multiblade  coordinate 
transformation  is  used  to  convert  the  rotor  states  from  the  rotating  to  the  fixed  coordinate 
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system.  The  controls  are  defined  in  a  rotating  coordinate  system.  The  state  matrix  A  and 
the  control  matrix  B  are  computed  using  numerical  perturbations  at  9  equidistant  azimuth 
steps  over  one  quarter  of  a  revolution,  corresponding  to  36  samples  over  one  revolution,  and 
are  expressed  in  the  form  of  a  Fourier  series.  For  the  time-periodic  case,  the  entire  Fourier 
series  is  used;  for  the  time-invariant  case,  only  the  constant  term  of  the  series  is  retained. 
Because  the  four  rotor  blades  are  assumed  to  be  identical,  the  Fourier  series  only  contains 
harmonics  that  are  integer  multiples  of  four.  In  this  chapter,  the  harmonics  retained  are  4/, 
8/  and  12/rev. 

3.6.1  Baseline  configuration 

Figures  3.1  shows  the  real  parts  of  the  poles  of  the  first  lag  mode  as  a  function  of  //,  computed 
using  the  constant  coefficient  approximation  (solid  lines)  and  the  rigorous  periodic  model 
(dashed  lines).  In  all  the  results  presented  in  this  section,  the  word  “pole”  means  “stability 
eigenvalue”  for  the  constant  coefficient,  or  LTI,  cases,  and  “characteristic  exponent”  for 
the  periodic,  or  LTP,  cases.  Also,  all  the  periodic  coefficient  results  are  obtained  using 
the  periodic  Schur  method.  The  figure  confirms  the  well  known  behavior  of  soft-in-plane 
rotor  stability,  namely,  that  the  low  frequency  lag  mode  is  the  least  damped.  This  mode 
is  traditionally  called  “regressive”  although,  strictly  speaking,  both  the  highest  and  lowest 
frequency  lag  modes  in  this  case  are  progressive  modes  (for  example,  see  [24]  for  further 
details).  Also,  hover  is  the  least  damped  flight  condition,  and  the  lag  damping  increases 
with  advance  ratio.  Finally,  as  expected,  for  the  case  of  very  small  values  of  the  advance 
ratio,  the  poles  computed  on  the  basis  of  Floquet  theory  coincide  with  the  poles  of  time- 
invariant  systems,  while  the  difference  increases  slightly  with  the  advance  ratio. 

Root  locus  plots  of  the  poles  and  zeros  of  the  first  lag  mode  with  advance  ratio  /i  as  the 
parameter  are  shown  in  Figures  3.2  and  3.3  for  the  LTI  and  LTP  cases  respectively.  It  should 
be  noted  that  the  zeros,  besides  being  dependent  on  the  specific  choice  of  input  and  output, 
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cannot  be  associated  with  normal  modes  in  the  same  sense  that  the  poles  can.  The  four  zeros 
in  the  figures  are  those  that  are  closest  to  the  four  poles  of  the  first  lag  mode,  and  are  given 
the  same  name  as  the  specific  pole  to  which  they  are  closest.  A  combination  of  pole  and 
zero  is  often  called  a  “dipole’' .  If  pole  and  zero  are  coincident,  then  the  dynamics  associated 
with  that  pole  are  canceled  by  the  zero,  and  therefore  they  do  not  appear  in  the  transfer 
function.  For  example,  in  Figure  3.2,  poles  and  zeros  of  the  progressive  and  regressive  lag 
modes  almost  exactly  cancel  one  another  in  hover.  As  a  consequence,  the  heave  acceleration 
response  to  collective  is  almost  unaffected  by  the  progressive  lag  mode,  and  it  would  be 
legitimate  to  remove  the  progressive  lag  mode  to  obtain  a  simplified  model  of  this  response. 
Similarly  for  the  regressive  lag  mode.  In  other  words,  pole-zero  cancellation  for  a  certain 
mode  decouples  that  mode  from  the  input-output  dynamics  being  considered.  Conversely, 
as  the  distance  between  pole  and  zero  of  a  dipole  increases,  the  corresponding  mode  plays 
an  increasing  role.  Therefore,  the  pole-zero  distance  can  be  seen  as  a  measure  of  coupling. 

Figure  3.2  shows  that,  in  general,  pole-zero  cancellation  does  not  occur,  and  that  the 
lag  dynamics  cannot  be  neglected  a  priori  in  computing  the  heave  acceleration  response  to 
collective  (the  extent  of  the  effect  of  lag  dynamics  will  depend  on  the  range  of  input-output 
frequencies).  The  figure  also  shows  that  the  zeros  of  one  of  the  collective  modes  and  of 
the  progressive  mode  become  less  damped  than  the  respective  poles,  although  they  remain 
minimum  phase  (i.e.,  with  negative  real  part).  Because  the  zeros  tend  to  “attract”  the  poles 
as  the  gains  of  a  feedback  loop  are  increased,  this  implies  that  the  closed  loop  system  may 
have  a  lower  lag  stability. 

Figure  3.3  shows  the  same  type  of  information  as  Fig.  3.2,  except  that  the  effect  of  peri¬ 
odic  coefficients  is  fully  taken  into  account,  and  therefore  the  “poles”  are  the  characteristic 
exponents  of  the  system.  The  same  general  considerations  as  for  Fig.  3.2  apply.  Compar¬ 
ing  Fig.  3.3  and  Fig.  3.2,  one  can  see  only  small  differences  between  the  behavior  of  the 
time-invariant  and  of  the  time-periodic  zeros.  This  indicates  that  the  constant  coefficient 
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approximation  is  generally  as  accurate  for  zeros  as  it  is  for  poles. 

Figure  3.4  compares  the  real  parts  of  the  LTI  and  LTP  zeros  nearest  to  the  lag  poles 
(each  zero  is  given  the  name  of  the  closest  pole).  As  is  the  case  for  the  poles,  in  hover 
the  LTI  and  LTP  zeros  essentially  coincide;  in  fact,  in  hover  the  only  periodicity  is  caused 
by  the  small  amount  of  cyclic  pitch  required  to  balance  the  tail  rotor  forces  and  moments. 
The  regressive  lag  zero  has  a  substantially  smaller  real  part  than  the  corresponding  pole, 
indicating  a  potential  reduction  of  inplane  damping  in  closed  loop  operation. 

Figure  3.5  shows  a  comparison  between  the  set  of  all  the  zeros  closest  to  the  imaginary 
axis  for  the  LTI  and  LTP  transfer  functions.  The  flight  speeds  range  from  hover  to  150  kts, 
corresponding  to  values  of  advance  ratio  from  0  to  0.35.  Similarly  to  the  time-periodic 
poles,  the  imaginary  parts  of  time-periodic  zeros  are  defined  within  a  multiple  of  the  period. 
Therefore  the  numerical  values  of  the  imaginary  parts  resulting  from  the  calculations  are  not 
necessarily  the  correct  ones.  Again,  similarly  to  the  time-periodic  poles,  each  time-periodic 
zero  can  be  identified  by  continuation,  i.e.,  by  starting  at  hover,  where  the  periodicity  is 
small,  and  time-invariant  and  time-periodic  poles  are  essentially  the  same,  and  following  the 
evolution  of  the  zero  as  /i  increases.  This  procedure  was  not  followed  for  the  time-periodic 
zeros  shown  in  Fig.  3.5,  and  therefore  the  imaginary  parts  of  those  zeros  are  not  necessarily 
meaningful.  On  the  other  hand,  the  real  parts  are  meaningful,  and  therefore  the  figure  shows 
the  presence  of  lowly  damped  or  nonminimum  phase  zeros.  Several  very  lowly  damped  and 
one  nonminimum  phase  zeros  are  clearly  visible. 

This  may  have  important  analysis  and  design  implications  from  a  control  engineering 
point  of  view.  In  general,  the  higher  the  frequency  of  a  stable  zero  near  the  imaginary 
axis,  the  higher  the  performance  limitations;  the  reverse  is  true  for  NMP  zeros,  i.e.,  the 
performance  limitations  are  higher  if  the  frequency  of  the  zero  is  lower. 
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3.6.2  High  thrust  configuration 


Figures  3.6  through  3.10  show  results  for  the  same  soft-in-plane  configuration,  but  at  a  higher 
thrust  condition,  i.e.,  with  a  value  of  ct/ct  =  0.1. 

Figure  3.6  shows  the  real  parts  of  the  poles  (stability  eigenvalues  for  the  time-invariant 
case,  and  characteristic  exponents  for  the  time-periodic  case)  of  the  first  lag  mode  as  a 
function  of  advance  ratio.  Compared  with  the  baseline  case,  Figure  3.1,  the  overall  damping 
of  the  mode  is  higher:  this  is  a  well  known  result,  which  depends  on  the  higher  aerodynamic 
damping  generated  by  the  higher  thrust. 

The  same  poles  and  the  four  zeros  closest  to  them  are  plotted  on  the  complex  plane  in 
Fig.  3.7  for  the  constant  coefficient  approximation,  and  as  a  function  of  speed.  From  the 
point  of  view  of  pole-zero  cancellation,  the  same  general  pattern  as  in  the  baseline  case, 
Fig.  3.7,  is  visible.  One  of  the  two  pole-zero  pairs  corresponding  to  the  collective  mode 
cancels  out  almost  exactly  at  all  speeds.  The  regressive  mode  pair  also  nearly  cancels  out. 
Significant  separations  can  be  seen  for  the  progressive  and  the  other  collective  mode  dipole. 
The  zero  of  the  collective  dipole  remains  very  close  to  the  imaginary  axis. 

The  same  type  of  information  is  shown  in  Fig.  3.8  for  the  time-periodic  case.  Comparing 
Figs.  3.7  and  3.8  shows  that  the  constant  coefficient  approximation  is  about  as  accurate  for 
the  zeros  as  it  is  for  the  poles.  The  only  exception  is  the  damping  of  the  progressive  lag 
zero  at  V  =  150  kts  and  /i  «  0.35,  which  is  substantially  underestimated  by  the  constant 
coefficient  approximation.  This  is  also  evident  from  Fig.  3.9,  which  shows  the  real  part  of 
the  same  zeros  as  a  function  of  the  advance  ratio  fi.  Finally,  Fig.  3.10  compares  LTP  and 
LTI  predictions  for  the  zeros  closest  to  the  imaginary  axis  and  in  the  right-half  plane.  As  in 
Fig.  3.5,  the  ambiguity  in  the  imaginary  parts  of  the  periodic  zeros  have  not  been  resolved, 
and  therefore  only  the  real  parts  are  always  correct.  The  same  general  arrangement  as  in  the 
baseline  case,  Fig.  3.5,  can  be  seen,  with  several  very  lowly  damped  and  one  nonminimum 
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phase  zero.  Therefore,  the  effect  of  thrust  on  the  character  of  the  zeros  does  not  appear  to 
be  significant. 

3.6.3  Stiff-in-plane  configuration 

Figures  3.11  and  3.12  refer  to  a  stiff-in-plane  configuration  obtained  by  modifying  the  lag 
stiffness  of  the  baseline  configuration  so  as  to  obtain  a  lag  frequency  in  hover  of  approximately 
1.4/rev.  Figure  3.11  shows  the  real  part  of  the  LTI  and  LTP  poles  of  the  lag  modes  as  a 
function  of  the  advance  ratio  /i.  As  typical  for  stiff-in-plane  hingeless  rotor  configurations, 
the  lag  degree  of  freedom  is  stable  in  hover,  but  its  stability  decreases  dramatically  at  high 
speed,  and  the  rotor  becomes  unstable  at  /x  ~  0.35.  There  are  some  small  differences  between 
LTI  and  LTP  predictions  at  hover:  these  are  caused  by  the  periodicity  introduced  by  the 
cyclic  pitch  required  to  balance  the  effects  of  the  tail  rotor. 

Figure  3.12  shows  the  real  parts  of  the  zeros  closest  to  the  four  lag  poles,  as  a  function 
of  fi.  The  corresponding  figures  for  the  soft-in-plane  case,  Figs.  3.4  and  3.9,  showed  four 
zeros,  whereas  Fig.  3.12  only  shows  three.  In  fact,  for  the  stiff-in-plane  case,  only  three 
zeros  are  so  close  to  corresponding  poles  that  they  can  be  clearly  identified  as  being  part 
of  a  dipole.  While  other  zeros  are  in  the  general  vicinity  of  the  fourth  lag  pole,  none  can 
be  associated  with  it  unambiguously.  Therefore,  because  there  is  not  even  an  approximate 
pole-zero  cancellation  for  this  lag  pole,  the  lag  mode  is  likely  to  participate  in  a  noticeable 
way  to  the  heave  acceleration  response  to  collective  pitch.  For  the  stiff-in-plane  case  too  the 
accuracy  of  constant  coefficient  approximation  begins  to  deteriorate  for  lower  values  of  /i  for 
the  zeros  than  for  the  poles.  This  is  shown  by  the  differences  between  the  exact  LTP  poles 
and  the  approximate  LTI  poles,  clearly  visible  in  Fig.  3.12,  which  become  noticeable  for  as 
low  as  /i  =  0.2.  Finally,  the  figure  shows  that  a  nonminimum  phase  zero  exists  at  hover. 
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3.7  Summary  and  Conclusions 


This  chapter  presented  a  methodology  for  the  calculation  of  the  zeros  of  transfer  functions  as¬ 
sociated  with  coupled  rotor-fuselage  systems,  considering  governing  equations  with  constant 
and  periodic  coefficients.  Only  the  single-input  single-output  case,  and  the  specific  transfer 
function  from  collective  pitch  to  CG  heave  accelerations  were  considered:  this  should  be  kept 
in  mind  when  generalizing  the  conclusions  of  this  study.  The  baseline  helicopter  configura¬ 
tion  was  similar  to  the  Eurocopter  BO-105.  A  flight  condition  with  higher  rotor  thrust  and 
a  stiff-in-plane  configuration  were  also  studied. 

The  main  conclusions  of  this  chapter  are: 

1.  The  zeros  of  a  coupled  rotor-fuselage  system  do  not  affect  its  open  loop  stability,  and 
therefore  can  be  ignored  if  no  active  control  system  is  present.  On  the  other  hand, 
if  an  active  control  system  is  present,  the  location  of  zeros  can  affect  the  closed  loop 
stability  characteristics,  and  may  trigger  closed-loop  instabilities  if  the  feedback  gain 
is  sufficiently  high.  The  presence  of  zeros  with  positive  real  parts  (nonminimum  phase 
zeros)  is  particularly  critical  in  determining  performance  limitations  of  the  control 
system. 

2.  Whereas  the  constant  coefficient  approximations  obtained  through  a  multiblade  coor¬ 
dinate  transformation  usually  give  accurate  predictions  of  the  stability  eigenvalues  for 
advance  ratios  of  up  to  n  —  0.3  —  0.35,  the  predictions  of  the  zeros  are  reasonably 
accurate  only  up  to  an  advance  ratio  of  about  /i  =  0.2,  and  occasionally  lower.  There¬ 
fore,  beyond  this  value  the  effects  of  periodic  coefficients  must  be  explicitly  taken  into 
account. 

3.  The  zeros  of  a  system  with  periodic  coefficients  can  be  computed  using  Floquet  the¬ 
ory,  however,  in  practice,  the  calculation  of  the  zeros  can  have  much  worse  numerical 
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properties  because  of  large  ratios  between  the  largest  and  the  smallest  zeros  (stiff¬ 
ness).  Therefore  the  conventional  shooting  techniques  normally  used  in  Floquet  sta¬ 
bility  analysis  often  prove  inadequate  and  more  accurate  computational  tools  must  be 
used.  These  large  zeros  cannot  simply  be  moved  to  infinity,  and  therefore  neglected, 
because  this  can  cause  inaccurate  predictions  of  the  zeros  near  the  imaginary  axis, 
which  in  turn  are  important  for  the  closed-loop  stability  of  the  coupled  rotor-fuselage 
system. 

4.  The  heave  acceleration  response  to  collective  pitch  is  affected  by  the  rotor  lag  dynamics. 
This  is  shown  by  the  fact  that  the  poles  of  the  lag  dynamics  are  not  canceled  by 
corresponding  zeros.  This  coupling  with  rotor  dynamics  is  stronger  at  higher  Ct/ct, 
and  even  more  for  stiff-in-plane  configurations. 

5.  The  transfer  function  from  collective  pitch  to  heave  acceleration  exhibits  several  zeros 
close  to  the  imaginary  axis  or  with  positive  real  parts.  The  number  and  location  of  these 
zeros  depend  on  the  configuration  and  flight  condition,  but  at  least  one  nonminimum 
phase  zeros  is  present  for  each  configuration,  and  at  every  speed.  These  zeros  can 
intrinsically  limit  the  benefits  achievable  from  active  control  systems,  and  therefore  it 
is  very  important  to  identify  them  as  part  of  the  design  of  any  such  control  system. 
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Figure  3.1:  Real  part  of  the  poles  of  the  first  lag  mode  as  a  function  of  /i;  baseline  configu¬ 
ration. 
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Figure  3.3:  LTP  lag  poles  and  zeros  as  a  function  of  /i. 
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Figure  3.4:  Real  part  of  LTI  and  LTP  lag  zeros  as  a  function  of  u. 


Figure  3.5:  Comparison  of  LTP  (left)  and  LTI  (right)  zeros  as  a  function  of  /i. 
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Figure  3.6:  Real  part  of  LTI  and  LTP  lag  poles  as  a  function  of  fi.  high  thrust  case 
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Figure  3.7:  LTI  lag  poles  and  zeros  as  a  function  of  /i,  high  thrust 
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Figure  3.8:  LTP  lag  poles  and  zeros  as  a  function  of  /i,  high  thrust  case. 
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Figure  3.10:  Comparison  of  LTP  (left)  and  LTI  (right)  zeros  as  a  function  of  /i,  high  thrnst 
case. 
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Figure  3.11:  Real  part  of  LTI  and  LTP  lag  poles  as  a  function  of  /i,  stiff-in-plane  case 
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Figure  3.12:  Real  part  of  LTI  and  LTP  lag  zeros  as  a  function  of  /i,  stiff-in-plane  case. 
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Figure  3.13:  Comparison  of  LTP  (left)  and  LTI  (right)  zeros  as  a  function  of  /i,  stiff- in-plane 
case. 
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Chapter  4 


Closed-Loop  Aeromechanical  Stability  of  Hingeless  Rotor  Helicopters  with 
Higher  Harmonic  Control 


4.1  Introduction 

Higher  Harmonic  Control  (HHC)  and  Individual  Blade  Control  (IBC)  have  been  considered 
for  many  years  as  a  viable  approach  for  the  design  and  the  implementation  of  active  rotor 
control  laws  aiming  at  the  attenuation  of  helicopter  vibrations  (see,  e.g.,  the  recent  survey 
papers  [20,  37]).  The  main  idea  of  HHC  and  IBC  is  to  try  to  attenuate  N/rev  vibratory 
components  in  the  fuselage  accelerations  (N  being  the  number  of  rotor  blades)  or  in  the 
rotor  hub  loads  by  adding  suitably  phased  N/rev  and  other  components  to  the  rotor  controls, 
either  in  the  fixed  (HHC)  or  rotating  (IBC)  frame.  A  number  of  studies  have  been  carried 
out  to  determine  the  feasibility  of  active  vibration  control  both  from  the  theoretical  and  the 
experimental  point  of  view;  in  particular,  as  far  as  the  analysis  of  the  dynamic  behavior 
of  the  single-input  single-output  (SISO)  HHC  is  concerned,  a  fundamental  result  was  given 
in  [40]  where  a  continuous  time  analysis  of  HHC  was  carried  out  for  the  first  time  and  it 
was  shown  that,  to  first  approximation,  the  classical  T-matrix  HHC  algorithm  (see  [35]) 
can  be  written  as  a  linear  time  invariant  dynamic  compensator.  More  recently,  however, 
it  has  been  proposed  to  try  to  exploit  the  interharm.onic  coupling  due  to  the  periodicity  of 
rotor  dynamics  in  forward  flight  to  achieve  the  attenuation  of  N/rev  vibrations  by  means  of 
lower  frequency  inputs,  such  as,  e.g.,  2/rev  or  3/rev  for  a  4  bladed  rotor.  To  this  purpose,  a 
generalization  of  the  T-matrix  algorithm  has  been  proposed  in  the  literature  (see  [32]),  but 
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no  detailed  theoretical  analysis  of  that  approach  has  been  carried  out  so  far.  As  the  above 
mentioned  generalization  of  the  T-matrix  algorithm  turns  out  to  be  a  linear  time-periodic 
compensator,  we  will  refer  to  it  as  the  Periodic  HHC  (PHHC)  algorithm.  Therefore,  both 
the  HHC  and  the  PHHC  algorithms  call  for  the  use  of  periodic  systems  theory  ([1])  for 
closed  loop  stability  and  performance  analysis.  However,  a  very  limited  attention  has  been 
devoted  so  far  in  the  literature  to  the  dynamic  analysis  of  vibration  attenuation  schemes; 
in  particular,  the  existing  contributions  to  the  study  of  closed  loop  stability  issues  (see  for 
example  [17])  deal  only  with  time-invariant  dynamic  models  of  helicopter  dynamics,  and  the 
assessment  of  the  role  of  periodicity  in  determining  the  actual  closed  loop  dynamics  still  has 
to  be  fully  assessed.  More  recently,  Cheng  et  al.  [9]  have  presented  a  methodology  for  the 
derivation  of  approximate  linearized,  time-invariant,  state-space  models  of  helicopters  and 
have  examined  the  interaction  between  HHC  and  FCS. 

In  the  light  of  the  above  remarks,  the  objectives  of  this  chapter  are  the  following: 

1.  to  provide  a  simple  state  space  derivation  for  the  continuous  time  form  of  the  SISO 
HHC  compensator  (input  and  output  at  the  same  rotor  harmonic),  first  introduced  in 
[40]; 

2.  to  demonstrate  how  the  same  approach  can  be  used  to  work  out  a  state  space  rep¬ 
resentation  for  the  SISO  PHHC  compensator  (input  and  output  at  different  rotor 
harmonics),  which  is  suitable  for  closed-loop  stability  and  robustness  analysis; 

3.  to  generalize  the  above  results  to  get  to  a  general  approach  for  the  derivation  of  the 
state  space  form  for  a  MIMO  HHC  controller  (any  combinations  of  harmonics  for  inputs 
and  outputs);  and 

4.  to  present  the  results  of  a  numerical  investigation  into  the  closed-loop  stability  prop¬ 
erties  of  Higher  Harmonic  Control,  based  on  a  simulation  study  of  the  coupled  rotor- 
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fuselage  dynamics  of  a  four  bladed  hingeless  rotor  helicopter. 


(Note  that  some  preliminary  results  have  been  presented  in  [27]). 

The  main  significance  of  this  work  is  that  it  addresses  the  question  of  whether  the  presence 
of  a  closed-loop  ffffC  will  affect  the  aeroelastic  stability  of  a  coupled  rotor-fuselage  system. 
Although  some  partial  answers  could  obviously  be  extracted  from  the  analysis  of  transient 
time  histories,  no  suitable  methodology  to  answer  this  question  directly  is  available  in  the 
published  literature. 

4.2  Helicopter  simulation  model 

The  baseline  simulation  model  used  in  this  chapter  is  a  nonreal-time,  blade  element  type, 
coupled  rotor-fuselage  simulation  model  (see  [38]  for  details).  The  fuselage  is  assumed  to  be 
rigid  and  dynamically  coupled  with  the  rotor.  A  total  of  nine  states  describe  fuselage  motion 
through  the  nonlinear  Euler  equations.  Fuselage  and  blade  aerodynamics  are  described 
through  tables  of  aerodynamic  coefficients,  and  no  small  angle  assumption  is  required.  A 
coupled  flap-lag-torsion  elastic  rotor  model  is  used.  Blades  are  modeled  as  Bernoulli-Euler 
beams.  The  rotor  is  discretized  using  finite  elements,  with  a  modal  coordinate  transformation 
to  reduce  the  number  of  degrees  of  freedom.  The  clastic  deflections  are  not  required  to  be 
small.  Blade  element  theory  is  used  to  obtain  the  aerodynamic  characteristics  on  each  blade 
section.  Quasi-steady  aerodynamics  is  used,  with  a  3-state  dynamic  inflow  model.  Linearized 
models  are  extracted  numerically,  by  perturbing  rotor,  fuselage,  and  inflow  states  about  a 
trimmed  equilibrium  position.  Because  the  equations  of  the  coupled  rotor /fuselage  dynamics 
are  written  in  the  fixed  frame  of  reference,  the  linearized  models  turn  out  to  be  time-periodic 
with  period  T /TV,  where  N  is  the  number  of  blades  and  T  is  the  period  of  one  rotor  revolution. 
Note  that  in  the  following  the  azimuth  angle  ijj  —  Of  will  be  used  as  independent  variable. 

The  matrices  of  the  linearized  model  are  generated  as  Fourier  series.  For  example,  the 
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state  matrix  A(ij})  is  given  as: 


K 

A{ip)  =  A0  +  ^2  [Akccos(kNifj)  +  Aks  sm(kNi/;)]  (4.1) 

k=\ 

where  the  matrices  A0,  Akc,  and  Aks  are  constant,  and  only  A0  is  retained  for  constant 
coefficient  approximations. 

Similarly,  the  control  matrix  B(gp)  is  obtained  assuming  for  the  pitch  control  of  each 
blade  the  form 

I< 

k=  1 

where  i  is  the  blade  number,  and  K  is  the  total  number  of  input  harmonics.  The  rotor  input 
vector  u  used  in  this  chapter  is 

U  =  [00  Ole  9\s  03c  Ozs  04c  04S  05c  05s]  (4-3) 

therefore,  for  the  4-bladed  rotor  considered  in  this  chapter,  the  higher  harmonic  input  is 
composed  of  the  N-l/rev,  N/rev,  and  N+l/rev  harmonics.  For  simplicity,  the  tail  rotor 
collective  input  90t  is  omitted  from  all  the  equations  of  this  chapter,  but  it  is  obviously 
included  in  the  actual  mathematical  model. 

As  Eq.  (4.2)  shows,  the  higher  harmonic  control  is  applied  in  the  rotating  system.  There¬ 
fore,  although  the  control  is  identical  for  all  blades  at  the  same  azimuth  angle,  this  arrange¬ 
ment  is  what  is  often  defined  as  Individual  Blade  Control  (IBC).  However,  it  should  be 
pointed  out  that  while  HHC  and  IBC  represent  significantly  different  technologies  from  the 
implementation  point  of  view  (e.g.,  choice  of  actuators  and  sensors),  they  are  completely 
equivalent  from  the  control  theoretic  point  of  view. 

Finally,  the  trim  procedure  is  the  same  as  in  Ref.  [8].  The  rotor  equations  of  motion  are 
transformed  into  a  system  of  nonlinear  algebraic  equations  using  a  Galerkin  method.  The 
algebraic  equations  enforcing  force  and  moment  equilibrium,  the  Euler  kinematic  equations, 


2,71 

9kc  cos  (  ip  T  —  i  )  +  9ks  sm 


,  2?r  • 
i,+  N‘ 


i  =  0, . . .  ,  N  —  1  (4.2) 
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the  inflow  equations  and  the  rotor  equations  are  combined  in  a  single  coupled  system.  The 
solution  yields  the  harmonics  of  a  Fourier  expansion  of  the  rotor  degrees  of  freedom,  the 
pitch  control  settings,  trim  attitudes  and  rates  of  the  entire  helicopter,  and  main  and  tail 
rotor  inflow.  During  the  trim  calculations  the  HHC  system  is  turned  off. 

4.3  State  space  formulation  of  higher  harmonic  controllers 

This  section  presents  the  derivation  of  the  state  space  formulation  of  a  MIMO  HHC  controller 
in  which  inputs  and  outputs  can  be  at  arbitrary  multiples  of  the  fundamental  rotor  frequency. 
This  is  done  by  building  state  space  formulations  of  HHC  controllers  of  increasing  complexity. 
After  some  background  on  the  T-matrix  algorithm,  a  continuous-time,  state  space  analysis 
is  presented  for  the  case  of  a  SISO  HHC  system  in  which  input  and  output  are  at  the  same 
harmonic  (in  this  case,  N/rev).  Next,  the  analysis  is  extended  to  the  case  in  which  input 
and  output  are  at  different  harmonics.  Finally,  the  case  of  the  MIMO  HHC  system  with 
inputs  and  outputs  at  arbitrary  harmonics  is  considered,  by  combining  the  results  of  the 
two  previous  cases.  More  precisely  the  following  three  cases,  corresponding  to  three  different 
selections  for  the  control  input  vector  u,  will  be  considered: 

•  Control  input  given  by  a  single  harmonic  at  the  blade  passing  frequency  N/rev: 

u  =  uN  =  [0Nc  9Ns]T 

•  Control  input  given  by  a  single  harmonic  at  a  frequency  M/rev  different  from  N/rev: 

u  =  uM  =  [ 9Mc  9Ms\ 

•  Control  input  given  by  the  superposition  of  a  number  of  different  harmonics,  such  as, 
for  example: 

u  =  [6/jv_i)c  9(n-i)s  9nc  9ns  9(n+ i)c  $(tv+i)s]  (4.4) 
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For  the  purpose  of  the  present  chapter,  only  fixed  parameters  HHC  will  be  considered, 
i.e.,  the  T-matrix  will  be  fixed  and  not  adaptive.  Also,  the  stability  analysis  will  be  carried 
out  in  continuous  time:  the  role  of  digital  implementation  on  the  stability  and  performance 
of  the  HHC  loops  will  be  investigated  in  future  work  [29]. 

4.3.1  Basic  T-matrix  algorithm 


A  typical  non-adaptive  HHC  system  is  based  on  a  discrete  time  mathematical  model  de¬ 
scribing  the  response  of  the  helicopter  to  higher  harmonic  inputs,  of  the  general  form 


yN{k)  =  TN:NuN(k )  +  ym(k) 


(4.5) 


where  k  is  the  rotor  revolution  index,  y at  is  a  vector  of  N/rev  harmonics  of  measured  outputs 
(e.g.,  hub  loads  or  accelerations  at  some  point  of  the  fuselage),  u at  is  a  vector  of  control  inputs, 
and  Tv,v  is  a  2  by  2  constant  matrix.  The  vector  y ^{k)  is  defined  as 


y  N(k)  = 


y  Nc(k) 
y  Ns(k) 


n  (  Aj— |—  1 )  7T 


^  J  kn 


y(fi>)  cos(N'ijj)  dfi 


1 
7 T 


^(fc+l)7T 


'  /C7T 


y (VO  sm(Nip)  dfi> 


(4.6) 


The  vector  yoAr  contains  the  N/rev  harmonics  of  the  “baseline”  vibrations,  i.e.,  the  vibrations 
in  the  absence  of  HHC.  In  practical  applications,  yw  is  typically  the  output  of  a  digital  or 
analog  harmonic  analyzer.  The  control  input  vector  is  similarly  defined  as: 


Uat 


9  Nc 

Ons 


(4.7) 


where  9nc  and  9ns  are,  respectively,  the  cosine  and  sine  components  of  the  N/rev  pitch 
control  input,  applied  in  the  rotating  system. 

The  HHC  inputs  are  generally  updated  at  discrete  time  intervals,  for  example,  once  per 
rotor  revolution.  The  conventional  HHC  control  law  is  derived  by  minimizing  at  each  discrete 
time  step  k  the  cost  function 


J  (k)  =  yN(k )J  QyN{k)  +  Aua^A^TAuaK*;)  (4.8) 
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where  Q  =  Q1  >  0,  R  >  0  and  Au a r(k)  is  the  increment  of  the  control  variable  at  time  k, 
i.e. 


AuN(k)  —  uN(k)  —  uN(k  —  1)  (4.9) 

Differentiating  Eq.  (4.8)  with  respect  to  Au N(k)  yields  the  control  law 

u N{k  +  1)  =  u N(k)  -  KN}NyN(k)  (4.10) 

where  KN>N  =  {Tj,  NQTN)N  +  i?)_1T^  NQ .  Equation  (4.10)  is  well  known  in  the  literature 
as  the  ’’T-matrix”  algorithm.  It  can  be  seen  from  Eqs.  (4.5)  and  (4.10)  that  this  control 
algorithm  introduces  a  discrete  time  integral  action  which  ensures  that  — »  0  as  k  — >  oo. 
Actually,  with  Q  =  /2, 2  and  R  =  0  deadbeat  control  (i.e.,  the  output  goes  to  zero  after  one 
discrete-time  step)  could  in  principle  be  achieved  if  exact  knowledge  of  the  TNjN  matrix  was 
available,  and  if  the  static  model,  Eq.  (4.5),  was  an  accurate  representation  of  rotor  dynamics. 
However,  these  two  assumptions  are  generally  not  satisfied,  as  TNtN  can  only  be  estimated 
and  Eq.  (4.5)  clearly  does  not  hold  if  the  helicopter  is  not  operating  in  steady  state.  Note, 
also,  that  if  in  the  cost  function,  Eq.  (4.8),  one  chooses  R  =  0  and  Q  proportional  to  the 
identity  matrix,  then  the  control  law,  Eq.  (4.10)  reduces  to 

uN(k  +  l)  =  uN(k) -T~'NyN(k)  (4.11) 

which  can  be  given  a  minimum  variance  interpretation,  in  the  sense  that  this  control  law 
guarantees  at  each  time  step  the  closed  loop  minimization  of  the  cost  function 

J(k)  =yN(k)TyN(k)  (4.12) 

Neglecting  the  effects  of  the  sample  and  hold  scheme  of  the  digital  implementation  in  the 
T-matrix  algorithm,  the  overall  control  algorithm  can  be  represented  by  the  block  diagram 
given  in  Figure  4.1. 
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Figure  4.1:  Block  diagram  of  the  continuous  time  SISO  HHC  algorithm. 

4.3.2  SISO  HHC  with  input  and  output  at  the  same  frequency 

Following  Ref.  [40],  choose  ]Jnc  and  y^s  as  state  variables  for  the  controller  in  Fig.  4.1.  Then, 
the  following  state  space  model  for  the  HHC  compensator  is  obtained: 


where 


y  n  =  AcyN  +  Bc(ip)  y 
u  =  Ccy  N 


Ac 

BcW 

Cc 


0  0 
0  0 


Kn,n 


cos(N'ip) 

sin(AT^) 


2 

T 


1-2,2 


4.3.3  SISO  HHC  with  input  and  output  at  different  frequencies 


(4.13) 

(4.14) 


(4.15) 

(4.16) 


The  HHC  input  in  the  rotating  system  is  usually  not  limited  to  the  same  N/rev  frequency  of 
the  vibrations  to  be  attenuated.  Typically,  inputs  at  N-l/rev  and  N+l/rev  are  also  applied 
(recall  that  N/rev  inputs  of  collective,  longitudinal,  and  lateral  cyclic  pitch  in  the  fixed 
system  result  in  N-l,  N,  and  N+l/rev  pitch  inputs  in  the  rotating  system). 
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In  this  case,  the  steady  state  model  relating  the  N/rev  harmonic  of  the  output  y(t)  to 
the  M/rev  harmonic  of  the  pitch  input  u(t),  with  M  ^  N,  can  be  written  in  the  form 

y  jv(fc)  =  TN,MuM(k)  +  y0N(k)  (4.17) 

where  u m  is  defined  as  in  Eq.  (4.7),  but  for  an  M/rev  harmonic,  and  where  the  (constant) 
matrix  Tm,m  relates  the  amplitude  of  the  M/rev  control  input  u  to  the  corresponding  steady 
state  amplitude  of  the  N/rev  component  of  the  output  y.  The  control  scheme  for  the  atten¬ 
uation  of  N/rev  vibrations  using  an  M/rev  input  can  then  be  derived  along  the  lines  of  the 
previous  case,  and  is  represented  by  the  equation 


u M{k  +  1)  =  u M(k)  -  KNjMyN(k)  (4.18) 

where  Kn,m  =  {Tn  mQTn,m  +  R)  1  TJjmQ.  As  shown  in  the  following  section,  the  matrix 
Tn,m  can  be  related  to  the  Harmonic  Transfer  Function  (HTF)  of  the  controlled  system, 
which  is  an  extension  to  periodic  systems  of  the  frequency  response  function  of  a  time- 
invariant  system  [41,  14].  In  addition,  as  in  the  case  of  HHC  with  input  and  output  at 
the  same  frequency  N/rev,  the  discrete  control  law,  Eq.  (4.18)  guarantees  that  — >  0  as 

k  — >  oo,  provided  that  the  system  can  be  modeled  as  in  Eq.  (4.17). 

Similarly  to  the  M  —  N  case,  the  state  space  model  for  the  case  N  ^  M  is  given  by 


where 


y  N  =  AcyN  +  Bc{ip)  y 
U  =  Ccy  N 
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(4.19) 

(4.20) 

(4.21) 

(4.22) 

(4.23) 
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This  discussion  shows  that  a  coupled  rotor- fuselage  system  with  even  a  simple  SISO 
HHC  controller  is  intrinsically  a  system  with  periodic  coefficients  if  the  HHC  output  and 
the  vibration  to  be  attenuated  are  at  two  different  multiples  of  the  rotor  frequency.  This 
happens  even  if  the  rotor-fuselage  system  is  modeled  as  a  system  with  constant  coefficients. 
Therefore,  rigorous  stability,  performance  and  robustness  analyses  of  an  HHC  system  can 
only  be  carried  out  using  the  tools  of  periodic  systems  theory. 

4.3.4  MIMO  with  input  and  output  at  arbitrary  harmonics 

In  typical  implementations  of  HHC,  multi-harmonic  signals  are  frequently  used  to  attenuate 
several  components  of  the  vibratory  loads.  For  example,  inputs  at  N-l/,  N/  and  N+l/rev, 
sine  and  cosine  (for  a  total  of  6  inputs),  could  be  used  simultaneously  to  control  six  compo¬ 
nents  of  the  N/rev  vibratory  hub  forces  and  moments.  Therefore,  this  section  extends  the 
previous  SISO  discussion  to  a  MIMO  HHC  system.  We  will  consider  a  general  configuration 
in  which  output  measurements  of  N/rev  vibration  are  available  at  n  different  locations,  while 
a  number  m  of  harmonics  at  frequencies  Ni,  i  =  1 , ,m  is  applied  on  the  control  input  u. 
In  this  case,  the  measurement  vector  has  2 n  elements  and  is  defined  as: 

Yn  =  [vnc  ■  ■  ■  Vnc  ■  ■  ■  Vns  ■  ■  ■  Vns]  (4.24) 

where  y) Vc  and  ylNs,  i  =  1, . . .  ,  n  are,  respectively,  the  cosine  and  sine  components  of  the  i-th 
N/rev  output,  which  can  be,  for  example,  a  force  or  moment  component,  or  a  component  of 
a  linear  or  angular  acceleration  at  one  or  more  points  of  the  fuselage. 

The  HHC  input  vector  u  has  2m  elements  and  is  defined  as 

UT  =  [uniC  UNlS  ■  ■  ■  Un2c  Un2s  ■  ■  ■  UNmC  UN^s)  (4-25) 

where  UNiC,  i  —  1,  ■  ■  ■  ,m  and  UNiS,  i  —  1,  ■  •  •  ,  m  are  the  cosine  and  sine  component  of  the 
HHC  input,  at  desired  harmonics  not  necessarily  equal  to  N/rev. 
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Assume  now,  as  in  the  SISO  case,  that  input  and  output  harmonics  are  related  by  the 


linear  equation 


y  N(k)  =  Tu(k)  +  y0N(k) 


(4.26) 


where  T  is  a  2 n  x  2m  constant  coefficient  matrix,  which  is  again  related  to  the  HTF  of  the 
time  periodic  linearized  model  of  the  helicopter.  Then,  the  ”T-matrix  algorithm”  is  given 
by 

u  (k  +  1)  =  u(/c)  —  Kyiy(k)  (4.27) 

where  K  =  ( TtQT  +  R)~1TtQ,  where  Q  =  QT  >  0  and  R  =  RT  >  0  are  cost  weighting 
matrices  of  suitable  dimensions. 

In  the  MIMO  case,  the  operation  of  the  HHC  control  law  differs  considerably  depending 
on  the  relationship  between  the  number  of  control  inputs  and  measured  variables  which  are 
available.  To  illustrate  this  difference,  the  formulation  of  the  ”T-matrix  algorithm”  in  the 
MIMO  case  with  Q  =  /n>n  and  R  —  0,  will  be  presented  separately  for  the  cases  of  n  =  m 
and  of  n  >  m.  (The  case  n  <  m,  i.e.,  with  more  controls  than  vibration  outputs  to  be 
attenuated,  is  not  likely  to  be  important  in  practice.) 

In  the  case  of  a  “square”  control  problem,  i.e.,  when  n  =  m,  the  SISO  algorithm  can  be 
readily  extended  to 

u  (k  +  1)  =  u  (k)  —  T~1yjv(k)  (4.28) 

On  the  other  hand,  if  n  >  m  the  matrix  T  is  no  longer  square  and  the  discrete  time  control 
algorithm  must  be  written  as 

u  (k  +  1)  =  u  (k)  —  T^yjy(k)  (4.29) 

where  T*  =  (TtT)~1Tt  is  the  pseudoinverse  of  T .  In  particular,  the  minimum  of  the  cost 
function  equals  zero  only  in  the  n  =  m  case,  i.e.,  unless  one  considers  (at  least)  the  square 
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case,  it  is  not  possible  to  guarantee  that  the  vibratory  disturbance  will  be  zeroed  on  all 
output  channels. 

The  equivalent  continuous  time  formulation  for  the  MIMO  HHC  compensator,  described 
in  discrete  form  by  Eq.  (4.27),  can  be  obtained  by  applying  the  previously  described  SISO 
results.  Therefore,  considering  first  the  case  of  a  control  system  with  as  many  inputs  as 
outputs,  the  state-space  formulation  is  given  by  the  order  2 n  system: 


Yn  =  AcYn  +  Bc(ip)  y 
u  =  Ccy  N 


where  Ar  is  the  2 n  x  2 n  matrix 


Ar  = 


0  0 
0  0 


0  0  ...  0 


Bc(ij))  is  the  2 n  x  n  matrix 


Bc{ij>)  =  K 


C0s(A^)4,r 
sin  (N^)Inir 


and 


Cc  —  r-pl2ny.2n 


(4.30) 

(4.31) 


(4.32) 


(4.33) 


(4.34) 


For  example,  consider  the  case  of  a  control  system  relying  on  the  application  of  N-l/,  N/, 
and  N+l/rev  inputs  in  the  rotating  frame  to  attenuate  one  component  of  the  accelerations 
in  n  —  3  different  locations  in  the  fuselage,  so  that  m  —  3,  N\  —  N  —  1,  N2  —  N  and 
^3  =  ^  +  1  and 


U1  —  [9(N-1)c  0(N-1)s  ^Nc  @Ns  9(N+1)c  ^(AT+1)s] 

yT  =  [lj\  V2  Va] 
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(4.35) 

(4.36) 


Then,  the  state  space  model  for  the  controller  is  given  by 


=  K 


(4.38) 


Ac  =  06x6  (4.37) 

[cos  (Ni/>)  0  0 

0  cos(AT/>)  0 

0  0  cos(Nip) 

sm(Nip)  0  0 

0  sin  (Nip)  0 
0  0  sin  (Nip) 

As  in  the  SISO  case,  because  the  control  inputs  are  directly  given  by  the  higher  harmonics 
of  6 ,  there  is  no  need  for  a  “modulation”  term  in  matrix  Cc  which  therefore  turns  out  to  be 
constant: 


Cc  —  ~7^hx6  (4.39) 

Similar  expressions  can  be  worked  out  in  the  case  of  a  control  system  with  more  outputs 
than  inputs. 


4.4  Definition  of  the  T  matrix  in  terms  of  the  helicopter  models 


The  control  laws  discussed  in  the  previous  section  call  for  the  availability  of  input/output 
models  for  the  helicopter  response  to  higher  harmonic  control  inputs.  This  section  provides 
some  background  on  the  frequency  response  of  time-periodic  systems  and  uses  such  analytical 
tools  to  derive  explicit  expressions  for  the  T  matrix. 

4.4.1  Development  of  the  Harmonic  Transfer  Function 

This  section  summarizes  the  main  aspects  of  the  development  of  the  Harmonic  Transfer 

Function  (HTF)  [41].  Consider  a  continuous-time  linear  periodic  system: 

x(t)  =  A(t)x.(t)  +  B(t)u(t) 
y  (t)  =  C(t)x(t)  +  D(t)u(t) 

Each  matrix  can  be  expanded  in  a  complex  Fourier  series 

OO 

A(t)=  Y,  4^ 

m=— oo 


(4.40) 

(4.41) 
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and  similarly  for  B(t ),  C(t)  and  D(t).  The  system  can  be  analyzed  in  the  frequency  domain 
as  follows.  Introduce  the  class  of  Exponentially  Modulated  Periodic  (EMP)  signals  [41].  The 
(complex)  signal  u(t)  is  said  to  be  EMP  of  period  T  and  modulation  s  if 

OO  OO 

u(t)  =  J2  u*eSkt  = eSt  UkejkQt  (4-42) 

k=— oo  k= — oo 

where  t  >  0,  Sk  =  s  +  jkfl,  and  s  is  a  complex  scalar. 

The  class  of  EMP  signals  is  a  generalization  of  the  class  of  T-periodic  signals,  i.e.,  of 
signals  with  period  T\  in  fact,  an  EMP  signal  with  s  =  0  is  just  an  ordinary  time-periodic 
signal. 

In  much  the  same  way  as  a  time  invariant  system  subject  to  a  (complex)  exponential 
input  has  an  exponential  steady-state  response,  a  periodic  system  subject  to  an  EMP  input 
has  an  EMP  steady-state  response.  In  such  a  response,  all  signals  of  interest  (x,  x,  y )  can 
be  expanded  as  EMP  signals.  By  deriving  Fourier  expansions  for  A(t),  B(t),  C(t )  and  D(t), 
it  is  possible  to  prove  that  the  EMP  steady-state  response  of  the  system  can  be  expressed 
as  the  infinite  dimensional  matrix  equation  with  constant  elements  [41] 

sX  =  (A-N)X  +  BU 

y  =  cx  +  vu 

where  X,  U  and  y  are  doubly  infinite  vectors  formed  with  the  harmonics  of  x,  u  and  y 
respectively,  organized  in  the  following  fashion: 

XT  =  [••■  xT2  x7],  x'(  x^  •••]  (4.44) 

and  similarly  for  U  and  y .  A,  B ,  C  and  T>  are  doubly  infinite  Toeplitz  matrices  formed  with 
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the  harmonics  of  >!(•),  /?(•),  C(-)  and  D(-)  respectively  as  follows 


A  = 


Ao 

A-i 

A-2 

A- 3 

A- 4 

Al 

Ao 

A-! 

A-2 

CO 

a2 

A\ 

Ao 

A-! 

A-2 

CO 

A2 

A, 

Aq 

A-i 

a4  a3  A2  A\  a0 


(4.45) 


(and  similarly  for  B ,  C  and  T>),  where  the  submatrices  An  in  Eq.  (4.45)  are  the  coefficients 
of  the  Fourier  expansion  of  matrix  A(t),  given  in  Eq.  (4.41).  To  relate  these  coefficients 
to  those  of  the  Fourier  series  expansion  in  trigonometric  form,  Eq.  (4.1),  recall  that  the 
Fourier  series  expansion  of  a  scalar  function  can  be  rewritten  in  complex  exponential  form, 
he.,  a(t)  =  a0  +  Y^T=  i  (anc  cos  nut  +  ans  sin  nut )  =  Y^k=-oo  akAkuJt,  with  ak  =  (akc-jaks)/ 2, 
and  a_fc  =  ( akc  +  jaks) / 2,  k  =  1,2, ... .  Therefore,  the  coefficients  of  Eqs.  (4.41)  and  (4.1) 
are  related  by: 

Ak  3 j 4fcs) 

k  =  1,2,...  (4.46) 

A-k  =  \{Akc  +  jAks) 

with  A0  identical  in  both  Eq.  (4.41)  and  Eq.  (4.1).  Similar  relations  hold  for  the  harmonics 
of  B,  C,  and  D. 

Matrix  J\f  is  a  block  diagonal  complex-valued  matrix  given  by 


-2/  0  0  0  0 

0  -I  0  0  0 

J\f  =  blkdiag{jnfl/}  =  j£l  \  - ■  0  0  0  0  0 

0  0  0  /0 

0  0  0  0  2/ 


(4.47) 


where  /  is  the  identity  matrix,  of  size  equal  to  the  number  of  states. 

From  Eq.  (4.43),  one  can  define  the  HTF  G(s)  as  the  operator: 

g(s)  d=  C[sl  -  {A  -  M)YlB  +  V  (4.48) 
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which  relates  the  input  harmonics  and  the  output  harmonics  (contained  in  the  infinite  vectors 
U  and  y  respectively).  Equation  (4.48)  is  the  extension  to  the  case  of  periodic  systems  of 
the  corresponding  constant  coefficient  expression  for  the  transfer  function 

G(s)  =  C[sI-A]~lB  +  D  (4.49) 

In  particular,  if  s  =  0,  which,  in  the  helicopter  case,  corresponds  to  the  steady-state  response 
of  the  system  to  a  periodic  input  of  basic  frequency  N/rev,  the  appropriate  input/output 
operator  for  periodic  systems  becomes 

0(0)  =  C[N  -  AYyB  +  V  (4.50) 

4.4.2  Definition  of  the  T  matrix 

The  Tn>n,  T/v,m  and  T  matrices  used  in  the  formulation  of  the  HHC  and  PHHC  algorithms 
can  be  related  to  the  elements  of  the  HTF  of  the  linearized  helicopter  model,  as  follows. 

First  of  all,  note  that  according  to  the  definition  of  the  control  input  vector  u  which  has 
been  adopted,  the  rotor  will  be  subject  to  a  proper,  steady  state  higher  harmonic  control 
input  whenever  the  control  vector  u  is  constant.  This  implies  that  to  define  the  T  matrix 
for  the  helicopter  we  only  have  to  study  the  response  of  the  periodic  helicopter  models  to 
a  EMP  input  with  s  —  0,  i.e.,  we  only  have  to  compute  the  input/output  operator  0(0). 
For  example,  consider  the  linear  time-periodic  system,  Eq.  (4.40),  and  the  constant  input 
u(t)  =  u0.  The  vector  U  corresponding  to  u (t)  =  u0  is  given  by 

UT  =  [••■  OOuJ  0  0'"]  (4.51) 

and  the  steady  state  response  y  of  the  periodic  system  is  given  by 

y  =  0(O)W  (4.52) 
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which  can  be  equivalently  written  as 


2/- 27V 

■  ■  ■  G —27V,— 27V  G—2N,—N  G-27V,0  G<-27V,7V  G^2N,2N  '  '  ' 

0 

y-N 

G-N-2N  G-N-N  G-Nfl  G-N,N  G-N,2N 

0 

yo 

= 

ry  ry  ry  ry  ry 

t*0,-2N  ^0-N  ^0,0  ^0,AT  ^0,2  N 

Mo 

yN 

Gn,-2N  Gn,-N  Gn,0  GtV.TV  Gn,2N 

0 

y2N 

G2N-2N  G2N,—N  G2n)0  G  27V, TV  G2Ny2N 

0 

From  Eq.  (4.53)  we  have  that 


(4.53) 


y-N 

1 

o 

_yN  _ 

Cb 

y 

o 

(4.54) 


and  converting  the  N/rev  harmonics  of  the  output  from  exponential  to  trigonometric  form 
we  have  that 


yNc 

—  o 

RealfG/vo] 

JjNs_ 

—  Zi 

_Imag[Gjv,0]_ 

(4.55) 


(note  that  G_jv,o  and  GV,o  are  complex  conjugate)  so  that  the  T  matrix  is  given  by 


T 


2 


Real  \G  /v,o] 
Imag[Gjv,o] 


(4.56) 


4.4.3  Construction  of  the  T  matrix 

From  a  practical  point  of  view,  the  above  theoretical  analysis  of  the  frequency  response  of 
periodic  system,  and  the  corresponding  definitions  for  the  T-matrix  relating  selected  input- 
output  frequencies  only,  rely  on  the  use  of  infinite  dimensional  matrices.  When  it  comes  to 
the  numerical  construction  of  the  T-matrix,  however,  one  has  to  resort  to  finite  dimensional 
approximations  of  the  system  matrices  A,  B ,  C,  and  T>.  Consider,  for  example,  the  problem 
of  constructing  the  T-matrix,  as  defined  in  Eq.  (4.56)  for  a  periodic  system  of  the  form  of 
Eq.  (4.40)  with  n  outputs,  m  inputs  and  ns  states.  First  of  all,  one  chooses  the  dimension 
of  the  expansions  A,  B,  C ,  and  V  for  the  state  space  matrices  A,  B ,  C ,  and  D,  in  terms  of 
the  number  of  block  rows  one  wants  to  take  into  account  in  A.  For  example,  if  we  choose 
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to  include  a  number  %  =  5  of  blocks  in  each  row  of  the  expansion  of  the  system  matrices, 
then  A  has  dimension  nsriB  x  nsns  and  is  given  by 


(4.57) 


and  similarly  for  B ,  C,  and  T>.  Therefore,  the  HTF  is  given  by  the  2 img  x  mriB  matrix,  as 
follows 


Ao 

A-! 

A-2 

41-3 

A- 4 

A, 

Ao 

A-i 

A-2 

41-3 

A  = 

A2 

A\ 

Ao 

A-i 

A-2 

A3 

A2 

Ai 

Ao 

41-i 

_a4 

A3 

A  2 

A\ 

4lo  . 

V-2N 

y-N 

yQ  =  0(0)  U  = 
yN 

l)2N 


G-2N-2N  G-2N-N  G-2N,Q  G-2N,N  G-2N,2N 

G-N-2N  G-n,-N  G-n,  o  G-N,N  G^n,2N 

S~1  /~i  /'"Y  /'"Y 

^0,-2^  ^0,-AT  <^0,0  ^0,iV  <-*0,2iV 

GV,  -2JV  Gn,  -N  GNfi  GNjn  GN‘2N 
GW-27V  G2N-N  GWo  G2N,N  G2N,2N 


(4.58) 


The  blocks  G_ jv,o,  Gat,o  can  be  extracted  from  0(0)  as  the  submatrices  Qij(0),i  =  2 n  + 
1, . . .  ,  3n  and  j  =  2m  +  1, . . .  ,  3m,  and  £7^(0),  i  =  4n  +  1, . . .  ,  5 n  and  j  =  2m  +  1, . . .  ,  3m 
respectively.  Clearly,  the  choice  of  the  number  of  block  rows  ub  will  affect  the  accuracy  of  the 
numerical  construction  (see  Ref.  [42]  for  an  analysis  of  the  effect  of  truncation  in  the  study  of 
frequency  response  operators),  so  as  a  general  rule,  tib  should  be  chosen  sufficiently  large  to 
ensure  that  the  T-matrix  constructed  from  the  truncated  ffTF  gives  a  good  approximation 
to  the  actual  T-matrix. 


4.5  Formulation  of  the  coupled  helicopter/HHC  model 


The  compensator  will  be  designed  along  the  lines  of  Ref.  [41].  Denote  with  A(ip), 

C(ip),  and  D{rjj)  the  matrices  for  the  LTP  state  space  model  of  the  helicopter,  for  the 

selected  input/output  pair.  Similarly,  denote  with  Ac(ip),  B c('0),  Cc{ij})  the  compensator’s 

state  space  model.  Then  the  closed-loop  LTP  state  matrix  Ae{ij. ;)  is  given  by 

a  =  [  AW  B(ip)Cc(ip) 

’  lBc{ii>)C(iP)  Ac{1>)  +  BJ^DWCcW) 

The  closed-loop  stability  of  the  helicopter  with  HHC  is  then  given  by  the  characteristic 
exponents  of  2le  (?/>),  and  will  be  studied  as  a  function  of  the  design  parameters  Q  and  R. 


(4.59) 
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4.6  Results 


This  section  presents  closed-loop  stability  and  response  results  for  a  coupled  helicopter-HHC 
system.  The  stability  results  are  obtained  from  the  linearized  model  of  Eq.  (4.59)  or  its 
constant  coefficient  approximation.  The  closed-loop  response  results  are  obtained  from  the 
full  nonlinear  simulation  of  the  coupled  helicopter-HHC  system. 

The  helicopter  configuration  used  for  the  present  chapter  is  similar  to  the  Eurocopter 
B0-105,  with  a  thrust  coefficient  Ct/c r  =  0.071.  Three  blade  modes  are  used  in  the  modal 
coordinate  transformation,  namely,  the  fundamental  flap,  lag,  and  torsion  modes,  with  nat¬ 
ural  frequencies  of  1.12/rev,  0.7/rev,  and  3.4/rev,  respectively.  Because  the  aerodynamic 
model  consists  of  a  simple  linear  inflow  with  quasi-steady  aerodynamics,  vibratory  loads  and 
CG  accelerations,  and  consequently  also  HHC  inputs,  tend  to  be  underestimated.  Therefore, 
their  absolute  values  can  be  considered  representative  only  in  a  qualitative  sense.  However, 
the  overall  simulation  model  is  likely  reasonable  for  stability  studies,  and  for  a  general  as¬ 
sessment  of  the  design  and  closed-loop  analysis  methodology. 

For  all  the  vibratory  response  results,  the  helicopter  is  first  trimmed  in  steady,  straight 
flight  at  the  desired  velocity  with  the  HHC  system  turned  off.  Then,  the  nonlinear  simulation 
begins,  with  the  pilot  controls  held  fixed  at  their  trim  values  and  the  HHC  system  turned 
on  at  time  t  =  0.  The  Q  and  R  matrices  in  Eq.  (4.8)  have  been  defined  as  Q  =  T6  G  and 
R  =  r/6j6,  respectively,  where  /6j6  is  an  identity  matrix  with  six  rows  and  columns,  and  r  is 
a  parameter  that  varies  from  r=0  (no  restriction  on  control  effort)  to  r=1000. 

4.6.1  Results  for  V=80  kts 

Figure  4.2  shows  the  peak-to-peak  magnitude  of  the  4/rev  component  of  the  vertical  (i.e., 
along  the  z-body  axis)  acceleration  at  the  CG,  for  a  speed  V  =  80  kts,  corresponding 
to  n  =  0.19.  The  figure  shows  four  curves,  one  each  for  values  of  r=0,  10,  100,  and  1000, 
corresponding  to  increasing  restrictions  on  the  control  effort.  The  high-frequency  oscillations 
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visible  in  the  curves  of  this,  and  of  many  subsequent  figures,  are  largely  an  artifact  of 
the  numerical  procedure  used  to  extract  the  4/rev  magnitude  and  phase  from  the  time 
histories  of  the  accelerations.  Clearly,  the  HHC  system  is  very  effective,  and  reduces  the 
4/rev  vertical  acceleration  to  a  small  fraction  of  its  trim  value  in  just  a  few  seconds.  The 
vibration  attenuation  is  also  very  clear  for  the  CG  roll  acceleration  p  and  pitch  acceleration 
q:  the  magnitudes  of  the  4/rev  components  are  shown  in  Fig.  4.3.  Both  p  and  q  are  reduced 
to  5%  or  less  of  their  trim  values  in  no  more  than  6-7  seconds. 

Figure  4.4  shows  the  magnitude  of  the  3/,  4/,  and  5/rev  components  of  the  HHC  input  for 
the  cases  r  —  0  and  r  =  1000.  Figure  4.5  shows  the  corresponding  phase  angles.  Comparing 
the  two  sets  of  results,  it  can  be  seen  that  the  controls  reach  their  steady-state  values  much 
more  quickly  for  the  case  r  —  0  than  for  r  =  1000.  In  the  latter  case,  the  steady-state  values 
of  #3  and  #5  have  not  yet  been  reached  at  the  end  of  7  sec  of  simulation. 

It  is  interesting  to  note  that  the  action  of  the  HHC  system,  and  the  consequent  vibration 
reduction,  occurs  within  times  of  the  order  of  5-7  sec  or,  equivalently,  of  about  1  rad/sec. 
These  are  also  typical  time  scales  for  flight  control  systems,  and  also  overlap  typical  piloting 
frequencies.  Therefore,  the  results  previously  shown  indicate  the  possibility  of  interaction 
with  the  stability  and  control  characteristics  of  the  helicopter. 

It  is  also  interesting  to  consider  the  closed-loop  eigenvalues  of  the  system.  The  com¬ 
putation  of  the  closed- loop  state  matrix  Ae,  Eq.  (4.59),  was  carried  out  by  linearizing  the 
augmented  nonlinear  set  of  equations.  Figure  4.6  shows  a  root  locus  plot  of  just  the  controller 
eigenvalues  for  increasing  values  of  r,  using  a  constant  coefficient  approximation  to  Ae.  The 
system  displays  a  mildly  unstable  complex  conjugate  pair  at  80  knots,  but  there  is  no  trace 
of  the  instability  in  the  closed-loop  simulations  using  the  full  nonlinear  system,  previously 
shown.  The  instability  is  probably  due  to  the  errors  made  in  modeling  the  periodic  system 
with  a  constant  coefficient  approximation.  In  fact,  when  the  periodicity  is  fully  taken  into 
account,  the  instability  disappears.  This  can  be  seen  in  Fig.  4.7,  which  shows  the  real  parts 
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of  both  the  stability  eigenvalues  of  the  LTI  system  and  the  Floquet  characteristic  exponents 
of  the  LTP  system,  for  the  least  damped  modes.  None  of  the  characteristic  exponents,  which 
include  controller,  rotor,  and  rigid  body  modes,  becomes  unstable  for  any  of  the  values  of 
r  considered.  This  confirms  that,  whenever  the  HHC  input  includes  harmonics  that  are 
different  from  the  harmonic  that  one  is  trying  to  attenuate,  the  closed-loop  problem  is  in¬ 
trinsically  time-periodic.  Constant  coefficient  approximations  may  not  yield  correct  correct 
closed-loop  stability  results,  as  in  this  case,  even  at  lower  advance  ratios,  where  constant 
coefficient  approximations  give  acceptable  results  for  the  open-loop  system.  Apart  from  the 
heading  eigenvalue  at  the  origin,  the  modes  with  the  lowest  amount  of  damping  are  those  of 
the  controller. 

Finally,  the  position  of  the  eigenvalues  appears  to  be  linked  to  the  vibration  reduction 
performance.  In  general,  for  the  highest  control  effort  (tuning  parameter  r  =  0)  controller 
eigenvalues  tend  to  be  farther  away  from  the  origin,  and  as  r  increases  they  come  closer  to 
it. 

4.6.2  Results  for  V=140  and  170  kts 

Figures  4.8  and  4.9  show  the  4/rev  CG  acceleration  components  at  a  speed  of  V=140  kts, 
corresponding  to  an  advance  ratio  /i  =  0.33.  The  magnitude  of  the  vertical  acceleration  is 
shown  in  Fig.  4.8.  The  HHC  is  extremely  effective,  and  reduces  the  magnitude  of  the  4/rev 
accelerations  to  almost  zero  within  about  7  seconds.  Near-perfect  attenuation  of  the  roll 
acceleration  p  can  be  seen  in  Fig.  4.9.  The  figure  also  shows  that  the  pitch  acceleration  q  is 
also  very  well  attenuated  by  the  HHC  system. 

The  magnitudes  of  the  corresponding  values  of  the  3/,  4/,  and  5/rev  inputs  are  shown  in 
Fig.  4.10  for  the  cases  r  —  0  and  r  =  1000.  The  steady-state  values  of  each  control  are  reached 
in  about  7  seconds,  therefore  the  time  scale  of  action  of  the  controller  is  approximately  the 
same  as  in  the  80  kts  case.  Differently  from  the  80  kts  case,  the  control  time  histories  for 
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r  —  0  and  r  =  1000  are  very  similar. 

The  LTI,  closed-loop  eigenvalues  for  V=140  kts  are  shown  in  Fig.  4.11.  At  this  speed, 
all  the  eigenvalues  are  stable,  with  the  partial  exception  of  a  complex  controller  eigenvalue, 
which  is  unstable  but  extremely  close  to  the  origin.  The  real  parts  of  both  the  stability 
eigenvalues  of  the  LTI  system  and  the  Floquet  characteristic  exponents  of  the  LTP  system, 
for  the  least  damped  modes,  are  shown  in  Fig.  4.12.  All  the  characteristic  exponents  are 
stable.  As  in  the  V=80  kts  case,  the  modes  with  the  lowest  amount  of  damping  are  those  of 
the  controller. 

Finally,  Figure  4.13  shows  one  result  for  the  case  V=170  kts,  corresponding  to  n  —  0.4. 
Note  that  the  simulation  cannot  determine  a  trim  state  at  this  speed.  Therefore,  the  drag 
of  the  fuselage  was  arbitrarily  reduced  until  a  trimmed  state  could  be  achieved.  Figure  4.13 
shows  baseline  and  HHC-on  magnitudes  of  the  4/rev  component  of  the  vertical  acceleration. 
Again,  the  HHC  is  very  effective  at  attenuating  vibrations,  and  the  attenuation  occurs  on 
the  same  time  scales  as  for  the  speeds  previously  shown.  Additional  results  were  obtained 
for  this  speed,  but  are  not  presented  here  for  reasons  of  space.  However,  the  overall  trends 
are  the  same  as  seen  for  the  V=80  kts  and  V=140  kts  cases,  except  that  the  closed-loop  LTI 
system  is  stable. 

4.7  Conclusions 

The  chapter  presented  the  development  of  a  state  space  formulation  for  a  Multi-Input  Multi- 
Output  (MIMO)  Higher  Harmonic  Control  (HHC)  system.  The  development  started  with  a 
simple  state  space  derivation  for  the  continuous  time  form  of  a  Single-Input  Single-Output 
HHC  compensator  with  input  and  output  at  the  same  rotor  harmonic;  then  the  same  ap¬ 
proach  was  extended  to  the  case  of  different  harmonics  in  input  and  output,  which  resulted 
in  a  periodic  SISO  HHC  compensator;  finally,  that  result  was  generalized  for  the  derivation 
of  the  state  space  form  for  a  MIMO  HHC  controller.  The  chapter  also  presented  results  of  a 
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numerical  investigation  into  the  performance  and  stability  properties  of  a  closed  loop  HHC 
system,  implemented  in  the  rotating  system,  based  on  a  simulation  study  of  the  coupled 
rotor-fuselage  dynamics  of  a  four  bladed  hingeless  rotor  helicopter. 

The  results  of  the  simulation  study  indicate  that: 

1.  The  HHC  controller  is  very  effective  in  reducing  the  desired  components  of  the  4/rev 
accelerations  at  the  center  of  gravity.  Because  the  aerodynamic  model  used  leads  to 
underestimating  these  vibratory  components,  the  absolute  values  of  the  reduction  and 
of  the  control  inputs  may  not  be  fully  reliable.  However,  the  percentage  reductions 
indicated  by  the  simulations  are  in  excess  of  80-90%. 

2.  The  vibration  attenuation  occurs  within  5-7  seconds  after  the  HHC  system  is  turned 
on.  This  is  equivalent  to  a  frequency  of  around  1  rad/sec,  which  is  within  the  frequency 
band  in  which  flight  control  systems  and  human  pilots  tend  to  operate.  Therefore,  the 
interactions  and  potential  adverse  effects  on  the  stability  and  control  characteristics  of 
the  helicopter  should  be  explored. 

3.  The  HHC  problem  is  intrinsically  time-periodic  if  the  HHC  inputs  include  frequencies 
other  than  the  frequency  one  wishes  to  attenuate.  This  is  true  even  if  the  rest  of 
the  model  is  assumed  to  be  time- in  variant.  In  these  cases,  the  closed-loop  stability 
results  obtained  using  a  constant  coefficient  approximations  may  be  incorrect  even  at 
lower  values  of  the  advance  ratio  //,  where  constant  coefficient  approximation  of  the 
open-loop  dynamics  are  accurate. 
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Figure  4.2:  Peak-to-peak  4/rev  vertical  accelerations  at  the  helicopter  center  of  mass  for 
V=80  kts  (/i  =  0.188). 
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Figure  4.3:  Peak-to-peak  4/rev  roll  (top)  and  pitch  (bottom)  accelerations  at  the  helicopter 
center  of  mass  for  V=80  kts  (/i  =  0.188). 
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Figure  4.6:  Root-locus  of  controller  eigenvalues  of  LTI  closed-loop  system,  V=80  kts. 


80 


Real  part  of  the  characteristic  exponent  (rad/sec) 


Figure  4.7:  Real  parts  of  eigenvalues  (top)  and  characteristic  exponents  (bottom)  of  the  least 
damped  modes,  V=80  kts. 
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Figure  4.9:  Peak-to-peak  4/rev  roll  (top)  and  pitch  (bottom)  accelerations  at  the  helicopter 
center  of  mass  for  V=140  kts  (//  =  0.330). 
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Figure  4.11:  Root-locus  of  controller  eigenvalues  of  LTI  closed-loop  system,  V=140  kts. 
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Figure  4.12:  Real  parts  of  eigenvalues  (top)  and  characteristic  exponents  (bottom)  of  the 
least  damped  modes,  V=140  kts. 


Chapter  5 


Basic  Concepts  in  the  Treatment  of  Discrete  Systems  with  Periodic  Coefficients 

5.1  State-space  and  input-output  models  for  periodic  systems 

The  state  space  description  of  a  periodic  system  in  discrete  time  is 


x(t  +  1) 

=  A(t)x(t )  +  B(t)u(t ) 

(5.1) 

y(t) 

=  C{t)x{t)  +  D{t)u{t) 

(5.2) 

where  x(t)  is  the  state  vector,  of  size  n,  y(t )  is  the  output  vector,  of  size  p,  and  u(t)  is  the 
input  vector,  of  size  m.  The  matrices  A(-) ,  B (■) ,  C (■) ,  and  D(-)  are  periodic  of  period  T. 

The  stability  of  the  system  can  be  assessed  by  analyzing  the  so-called  m.onodromy  matrix. 
Precisely,  let  t)  (t  >  r)  be  the  system  transition  matrix,  i.e. 

A(t-l)A(t-2)  t  >  t  (53) 

The  monodromy  matrix  is  defined  as  the  transition  matrix  over  one  period,  e.g.,  [r,  r  +  T— 1], 
and  is  denoted  by  Ta(t)  =  $a{t+T,  t).  Its  eigenvalues  do  not  depend  upon  r  and  are  named 
characteristic  multipliers.  The  system  is  stable  if  and  only  if  the  characteristic  multipliers 
lie  in  the  open  unit  disk. 

In  the  same  way  as  for  time-invariant  systems,  the  external  properties  of  a  periodic  system 
can  also  be  studied  in  a  fully  input-output  context.  The  basic  causal  relationship  supplies 


the  output  y(t) as  a  linear  combination  of  past  values  of  the  input  up  to  time  r: 

y[t)  =  M0(t)u(t)  +  Mi(t)u[t  —  1)  +  M2{t)u(t  —  2)  +  M3(t)u(t  —  3)  +  . . . 

OO 

=  '52  ~  j)  (5-4) 

3= 0 

The  matrix  coefficients  =  0, 1,2, ,  are  periodic  functions  of  period  T,  known  as 

periodic  Markov  coefficients.  These  Markov  parameters  are  linked  to  the  impulse  response 
of  the  system  in  a  way  that  can  be  easily  assessed  by  making  reference  to  the  simple  case 
where  the  input  is  a  scalar  variable  (m  =  1).  Indeed,  denoting  by  5(t)  the  impulse  function, 
i.e., 

*<*>  =  {  l  otherwise  <5-5> 

and  by  y^p(t)  the  response  of  the  system  at  the  impulsive  input  u(t)  —  S(t  —  i)  it  follows 
from  Eq.  (5.4)  that 

Vimpit)  =  (5.6) 

In  the  general  multi-input  case,  the  j-th  column  of  the  Markov  coefficients  Mt_i(t )  represents 
the  output  response  of  the  system  to  an  impulse  applied  at  time  i  to  the  j-th  component  of 
the  unit  vector. 

The  periodicity  of  the  Markov  coefficients  entails  that  the  output  response  of  the  system 
at  a  generic  time  instant  t  =  kT  +  s,  where  s  is  an  integer  such  that  0  <  s  <  T  —  1,  can 
be  written  as  a  finite  sum  of  the  output  responses  of  T  time  invariant  systems  indexed  in 
the  integer  s.  As  a  matter  of  fact,  consider  again  Eq.  (5.4),  and  evaluate  y{t)  at  the  instant 
t  =  kT  +  s.  We  have 

OO  OO 

y(t)  =  ^2Mj(t)u(t  -  j)  -»•  y(kT  +  s)  =  ^2Mj(kT +  s)u(kT +  s- j)  (5.7) 

3=0  3=0 

Focus  now  on  the  summation  in  Eq.  (5.7).  The  summation  index  j  counts  the  sample  points 
as  time  progresses.  The  index  j  can  be  rewritten  in  terms  of  the  period  T  as  j  =  IT  +  i, 


89 


where  t  counts  the  number  of  periods,  and  i  is  the  sample  point  within  the  period.  For 
example,  assume  that  a  period  is  composed  of  T  =  10  time  points.  Then  j  =  22  can  be  seen 
as  either  the  22nd  time  point  overall,  or  the  second  point  of  the  third  period,  i.e.,  22  =  2 
(number  of  periods)  •  10  (length  of  one  period)  +  2  (2nd  point  in  the  period),  and  therefore 
£  =  2,  T  =  10,  and  i  =  2.  With  this  in  mind,  one  can  replace  j  with  £T  +  i  in  Eq.  (5.7),  and 
replace  the  summation  over  j  with  a  double  summation  over  the  number  of  points  in  one 
period  and  over  the  number  of  periods,  i.e., 

OO  T—  1  OO 

£  =  ££  (s-s) 

j= o  i=o  e=o 

Making  these  substitutions  in  Eq.  (5.7)  one  obtains 


OO 

y(kT  +  s)  =  Mj(kT  +  s)ujkT  +  s  —  j) 
j= o 

T—l  oo 

=  ££  Mi+er(kT  +  s )  u[kT  +  s  —  (i  +  £T)] 
i= o  e=o 

T—l  oo 

=  ££  Mi+er(s)  u[kT  +  s  -  (i  +  £T )] 
i= o  e=o 

because  M  is  periodic,  i ,e.,Mi+er(kT  +  s)  — 

T—l  oo 

5^  J2Mi+eT(s)u[kT  +  s-i~£T)] 
i= o  e=o 


That  is, 

T—l 

y(kT  +  s)  =  ^  Vi,s(k) 
i= 0 


(5.9) 


(5.10) 


y,Jk)  =  W  Mi+tT(s )  u\kT  +  s-i  -  ('!')} 

OO 

=  J2Mi+tr{s)u[{k-e)T  +  s-i)]  (5.11) 

£=0 
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Now,  if  we  define 


UitS(k)  =  u(kT  +  s  —  i)  (5-12) 

then 

UitS{k  —  £)  —  u[(k  —  t)T  +  s  —  i)]  (5.13) 

which  appears  in  Eq.  (5.11),  which  can  therefore  be  rewritten  as 

OO 

Vi,s{k)  =  ^2  Mi+tT{s)uiiS{k  -  £)  (5.14) 

1=0 

From  these  expressions,  it  is  apparent  that  yi,s(k)  is  the  output  of  a  time-invariant  system 
having  Mj(s),  M;+r(s),  Mi+2r(s), . . . ,  as  Markov  parameters.  Note  the  role  of  the  different 
indexes  appearing  in  these  expressions:  s  is  the  chosen  tag  time  index  for  the  output  variable, 
s  —  i  is  the  tag  time  index  for  the  input  variable,  and  k  is  the  sampled  time  current  variable. 
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Example:  Use  Eq.  (5.9)  with  k  —  2,  T  —  4,  s  —  1,  and  kT  +  s  =  9. 


y(9) 


T— 1  oo 


ee  Mi+er{s)  m[/;T  +  s  —  i  —  £T)\ 

i=0  ^=0 

T-l  oo 

EE  Mi+4*(1)  w(9  -  %  -  £T) 


i= 0  £=0 
3 


E 

i=0 


Mj(l)  u( 9  -  i)  +  Mi+4(  1)  w(9  -  i  -  4)  + 


^=o 


^=i 


Mi+8(1)  «(9  -  i  -  8)  +  Mm2(l)  «(9  -  i  -  12)  +  . 


e=2 


1=3 


Y  Mi(l)  u( 9  -  i)  +  Y  m*+4(1)  «(5  -  i)  + 


i=0 


f=0 


y,  Mj+8(1)  w(l  -  i)  +  y  Mj+i2(l)  w(-3  -  i)  + 

z=0  7=0 

M0(l)w(9)  +  Mi(1)m(8)  +  M2(1)m(7)  +  M3(l)w(6)  + 
1144(1)^(5)  +  M5(1)m(4)  +  M6(l)w(3)  +  _M7(l)ii(2)  + 


M%{X)u(V)  +  Afg(l),u(0)  +  Mio(l)w( — 1)  +  Mn(l)w( — 2)  + 

M\2iX)u{~ ’3)  +  IVjT  1 3 ( 1 ) (  4)  +  M14(l)w(— 5)  +  A/45(1)m(— 6)  +  . . . 


M0(l)w(9)  +  M4(1)m(8)  +  M2(1)m(7)  +  M3(l)w(6)  + 
M4(1)u(5)  +  M5(1)u(4)  +  M6(1)w(3)  +  _M7(l)it(2)  + 
1U8  (1)^(1)  +  lUg  ( 1 )  (  0  ) 


because  we  assume  that  the  input  starts  at  time  7  =  0,  and  therefore  u{t)  =0  for 
Compare  Eq.  (5.15)  with  the  results  of  using  directly  Eq.  (5.4): 

OO 

y(t)  =  Y  MAt)u(kt^  i) 

3= 0 

=  M0(l)u(9)  +  M4(1)w(8)  +  M2(l)u(7)  +  1143(1)-u(6)  +  Af4(l)u(5)  + 
M5(1)m(4)  +  M6(l)u(3)  +  M7(1)u(2)  +  1W8(1)m(1)  +  A/g(l)u(0) 


(5.15) 
t  <  0. 


(5.16) 
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The  two  results  are  identical  (this  concludes  the  example). 

For  each  i  —  0, 1,  2, . . .  ,  T  —  1,  (i.e.,  for  each  time  point  in  the  period  T  for  the  input), 
and  each  s  —  0, 1,  2, . . .  ,  T  —  1,  (i.e.,  for  each  time  point  in  the  period  T  for  the  output),  one 
can  define 


Hi(z,s)  =f  J2M^+i(s)z~e 


£=0 


—  Mi(s)  +  Mi+x{s)z  +  Mi+2r{s)z  +  . . 

—  -Wj(s)  +  o  ~ I-  •  •  • 


(5.17) 

(5.18) 

(5.19) 


2)  Z* 

This  is  the  transfer  function  from  the  input  u^s{k)  to  the  output  yi,s(k)  (both  signals  seen 
as  functions  of  k).  To  see  this,  start  by  rewriting  Eq.  (5.14) 


Vi,a(k)  =  y ^  Mi+er(s)ui,s(k  -  £) 


Eq.  (5.14)  repeated 


e=o 


=  Mi(s)uijS(k )  +  Mi+T(s)uitS(k  -  1)  +  Mi+2T(s)ui,s(k  -  2)  +  . . . 


(5.20) 


From  the  basic  properties  of  the  ^-transform  we  have 

Ui,s(k  - 1)  =  z~luiiS{k)  (5.21) 

Ui,s(k  -  2)  =  z~2ui)S{k)  (5.22) 


Substituting  into  Eq.  (5.20)  we  have 


ViAk ) 


Mi(s)uiAk)  +  Mi+T(s)uiAk  -  1)  +  Mi+2T{s)uiAk  -  2)  +  . . . 

Mi(s)uiAk )  +  Mi+T(s)z~1uiAk)  +  Mi+2T{s)z~2UiAk )  +  •  •  • 


-  Mi+T(s)z  1  +  Mi+2T(s)z  2  +  . . .  ]  ui:, 

- V* - y 


m 


- V - 

d=  Hi(z,s)  (see  Eq.  (5.18)) 


(5.23) 


which  confirms  that  Hj{z,s )  is  the  transfer  function  UiAk)  to  yijS(k).  By  using  z  as  the 
one-shift  ahead  operator  in  time  k  (namely,  the  T-shift  operator  in  time  t),  and  resorting  to 
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a  mixed  z/k  notation,  one  can  connect  the  outputs  y(t)  to  the  inputs  u(t)  directly,  without 
going  through  the  intermediate  quantities  y(t)  and  u{t).  To  see  this,  begin  by  rewriting 
Eq.  (5.10) 

T—l 

y(kT  +  s)  =  Eq.  (5.10)  repeated 

i=0 

expand  the  summation 

=  yo,sik)  +  yi,s{k)  +  y2,s{k)  +  . . .  +  yT-i,s(k) 

use  Eq.  (5.23) 

=  H0(z,  s)u0,s(k)  +  Hi(z,s)uitS(k)  + 

+H2{z,  s)u2,s(k)  +  . . .  +  Ht-^z,  s)uT_hs(k)  (5.24) 

Now  recall  that  it  was 


from  which 


UitS(k)  =  u[kT  +  s  —  i)  Eq.  (5.12)  repeated 

u0  >s(k)  =  u{kT  +  s) 
ui,s(k )  =  u{kT  +  s  —  1) 

U2,s(k )  =  u(kT  +  s  —  2) 


uT~i,s(k)  =  u(kT  +  s  —  (T  —  1)) 

Substitute  into  Eq.  (5.24)  to  get 

y(kT  +  s)  =  Hq(z,  s)u{kT  +  s)  +  Hi(z,  s)u(kT  +  s  —  1)  + 

+H2(z,  s)u(kT  +  s  —  2)  +  . . .  +  ,  s)u(kT  +  s  —  T  +  1)  (5.25) 

but  we  would  like  to  define  y(kT  +  s)  in  terms  of 

u{kT  +  s),  u(kT  +  s  +  1),  u(kT  +  s  +  2), . . ..  ,  u(kT  +  s  +  T  —  1) 
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so  we  can  then  write  input-output  relations  at  the  same  time  points,  therefore  we  need  to 
transform  the  terms  u(kT  +  s— i)  into  u(kT  +  s+i).  To  do  this,  start  by  writing  Eq.  (5.25) 
in  a  different  order 

Eq.  (5.25)  Eq.  (5.25)  rearranged 

y[kT  +  s)  =  H0(z,  s)u(kT  +  s)+ 

+Hi(z,  s)u{kT  +  s  —  1)  +Ht~i(z,  s)u{kT  +  s  —  T  +  1) 

+H2(z,  s)u{kT  +  s  —  2)  +HT-‘i{z1  s)u(kT  +  s  -T  +  2) 

2(2,  s)u(kT  +  s  —  T  +  2)  -\-H2(z,  s^uikT  +  s  —  2) 

+Ht-2(z,  s)u(kT  +  s  —  T  +  1)  +Hi{z,  s)u(kT  +  s  —  1) 

+H0(z,  s)u(kT  +  s) 

Then  use  the  one-shift  ahead  operator  in  time  z,  u(s  —  T)=  z~1u(s): 

u(kT  +  s-(T-  1))  ->  u(kT  +  s  +  l-T)  =  z~'u{kT  +  s  +  1) 

u(kT  +  s  —  (T  —  2))  — ■>  u{kT  +  s  +  2  —  T)  =  z~lu{kT  +  s  +  2) 

u{kT  +  s  —  (2))  — ■>  u(kT  +  s  —  2)  =  z~lu(kT  +  s  —  2  +  T) 

w(A:T  +  s  —  (1))  — >  u{kT  +  s  —  1)  =  z~lu{kT  +  s  —  1  +  T) 

and  rewrite  Eq.  (5.25)  rearranged 

y(kT  +  s)  =  Hq(z,  s)u(kT  +  s)  + 

+Ht-i(z,  s)z~lu[kT  +  s  +  1) 

-\-Hrp—2{z,  s)z  u(kT  +  s  +  2) 

-\-H-2(z,  s)z~1u(kT  +  s  +  T  —  2) 

+Hi(z,  s)z~1u(kT  +  s  +  T  —  1)  (5.26) 

The  transfer  function  H,(z,  s )  will  be  referred  to  as  the  sampled  transfer  function  at  tag  time 
s  with  input-output  delay  i. 

The  input-output  periodic  model  is  said  to  be  rational  if  all  transfer  functions  (5.17) 
are  indeed  rational,  i.e.,  they  are  transfer  functions  of  finite-dimensional  (time-invariant) 
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systems.  In  this  case,  the  periodic  system  can  be  given  a  state-space  finite-dimensional 
realization,  as  shown  in  Ref.  [13].  The  problem  of  finding  a  minimal  periodic  realization  in 
the  form  of  difference  equation  is  investigated  in  Refs.  [3]  and  [26]. 

Introduce  now  the  unit  delay  operator  cr(-),  which  moves  forward  a  signal  by  one  step, 

i.e., 

au(t)  =  u(t  +  1)  a~lu(t)  =  u(t  —  1)  a~ku(t)  —  u(t  —  k)  (5. 27) 


Using  cr(-),  Eq.  (5.4)  can  be  given  a  compact  form: 


y(t)  =  M0(t)u(t )  +  Mi(t)u(t  —  1)  +  M2(t)u(t  —  2)  +  M3(t)u(t  —  3)  +  . . .  Eq.  (5.4) 
=  M0(t)u(t )  +  Mi(f)a_1w(f)  +  M2(t)a~2u(t)  +  M3(t)a~3u(t)  +  . . . 

—  [M0(t)  +  1(-)  +  M2(t)a  2(-)  +  M3{t)<7  3(-)  +  . . .  J  u(t)  (5.28) 

=  G(a,-)\,u(t)  (5.29) 

where 


G(<J, -)|t  —  M0(t)  +  x(-)  +  M2(t)a  2(-)  +  M3(t)a  3(-)  +  . . . 

is  the  input-output  transfer  operator.  The  operator  G(a,-)\t  is  periodic,  i.e., 

G(a,  ■) \t+T  =  G(a,  •) |t  for  every  t 


in  fact, 


(5.30) 


(5.31) 


G{(J,-)\t+T  —  [Mo(t  +  T)  +  M\{t  +  T)a  1(-)  +  Af2(t  +  T)a  2(-)  +  ...] 
—  \_Mo(t)  +  M\{t)a  1  (■)  +  AI2(t)a  2(-)  +  ...] 


=  G(a,-)\t 

because  the  Markov  coefficients  Mft)  are  periodic.  Moreover,  G(a,-)  |t  enjoys  a  pseudo- 
commutative  property  with  respect  to  the  delays.  Precisely 

a~k  G(a,-)\t—  G(a,-)\t_ka~k  for  every  t  (5.32) 
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To  see  this,  write 


o-  k  [M0(t)u(t)\ 


a  k  Mi(t)a  1u(t) 


a  k  M2(t)a  2u{t) 


M0(t  —  k)u{t  —  k) 
M0(t  —  k)a~ku{t) 

a~k  [Mi(t)u(t  —  1)] 
Mi(t  —  k)u(t  —  1  —  k) 
M\[t  —  /c)cr_fc_1-u(t) 
Mi(t  —  k)a~l  a~k  u(t) 

a~k  [M2(t)u(t  -  2)] 
M2(t  —  k)u(t  —  2  —  k) 
M2(t  —  k)a~k~2u{t) 
M2(t  —  k)o~2  a~k  u(t) 


(5.33) 


(5.34) 


(5.35) 


Then 


a  k  G(a,  -)\tu(t)  =  a  k  [M0(t)  +  Mi(£)cr  l(-)  +  M2(t)a  2(-)  +  . . .  ]  u(t) 

=  [_M0(t  —  k)  +  M\{t  —  k)a~l  +  M2(t  —  k)a~2  +  . . .  ]  cr~ku(t ) 

=  G(a,-)\t_ka~ku(t ) 

which  implies 

a~k  G(a,  ■)  \t  —  G(a,  ■) \t_k  a~k 
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Because  of  this  property,  one  can  rewrite  Eq.  (5.30)  as 


G(a,-)\t  =  M0(t)  +  1(-)  +  M2(t)a  2(-)  +  M3(t)a  3(-)  +  . . .  Eq.  (5.30)  repeated 

=  M0(t)  +  u  1 Mi(t  +  1)(')  +  &  2M2(t  +  2)(-)  +  <7  3 M%(t  +  3)(-)  +  . . .  (5.36) 

Hence,  the  operators  cr~k  and  G(a,  •)  do  commute  if  and  only  if  the  integer  K  is  multiple  of 
the  period  T. 

If  the  transfer  operator  is  evaluated  in  a  specific  time-point,  say  t,  it  results  in  a  transfer 
function,  henceforth  denoted  by  G(a,t).  In  view  of  Eq.  (5.17),  such  transfer  function  can  be 
written  as 

oo  T- 1 

G(a,  t)  =  Mk(t)a~k  =  Y  Hi  (aT ,  t)  a~l  (5.37) 

k=0  i= 0 

In  fact, 

OO 

G{a,t)  =  J2M^a~k 

k= 0 

rewrite  k  as  jT  +  i 

T—  1  oo  T—  1  oo 

=  5Z  S  MjT+i(t)ajT+t  =  YY1  Mn+i(t)  {vT)  J  v~l 
i= 0  j= 0  i= 0  j= 0 

rewrite  t  as  kT  +  s 

T- 1  oo 

=  5Z  S  Mn+i(kT  + s )  (ffjr) 

*=0  j=0 

but  because  M  is  periodic,  M]T+l ( kT  +  s)  =  MjT+1(s),  and  therefore 

T- 1  oo 

=  X)  S  M^+*(S)  (<tT)  ^ 

i=0  j= 0 

finally,  use  the  definition,  Eq.  (5.17),  Hi(z,s )  =  YlJLo  ^]T+% (s)z ,  with  z  =  a1 

T-l 

=  (5.38) 

i=0 

which  is  the  right-hand-side  of  Eq.  (5.37).  Also  note  that  a  generates  a  shift  of  one  time 
step,  whereas  z  generates  a  shift  of  one  period.  Therefore,  if  the  period  is  composed  of  T 
time  steps,  to  make  one  shift  of  Z  it  takes  T  shifts  of  a,  i.e.,  z  =  aT . 
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We  now  introduce  the  symbol  0  which  will  be  often  used  in  the  chapter: 


0  =  exp 


+  %  sin 


(5.39) 


Hence  1,  0,  02, . . .  ,  0T  1  are  the  T-roots  of  the  unit.  Note  that 

—  Y^0sfc={  1  s  =  0,±T,  ±2T, ... 

T  '  1  0  otherwise 

k= o  k 


(5.40) 


Then,  it  is  easy  to  see  that  one  can  recover  the  sampled  transfer  functions  Hi  (aT,  t)  from  a 
single  transfer  function  G(a,t)  as  follows: 


Hi(aT,t) 


1 

T 


T-l 


^  G(a0fc,  t)(pki 


,k= 0 


a 


(5.41) 


5.2  Time-lifted  reformulation 


In  this  section  we  introduce  the  most  classical  time-invariant  reformulation,  namely  the 
lifted  reformulation.  The  simplest  way  to  explain  the  idea  of  lifting  is  to  make  reference  to  a 
(discrete-time)  signal,  say  v(-).  Associated  with  v(-)  one  can  introduce  the  augmented  signal 


vr(k) 


v(kT  +  r) 
v{kT  +  r  +  1) 
v(kT  +  r  +  2) 

v{kT  +  r  +  T  —  2) 
v{kT  +  t  +  T  —  1) 


(5.42) 


where  r  is  a  tag  point,  arbitrarily  chosen.  We  will  also  use  the  symbol  v^\k)  for  the 
components  of  vector  vT(k),  i.e., 


v^\k)  =  v{kT  +  r  +  i  —  1)  %  —  1, 2, . . .  ,  T  (5.43) 


Note  that  v^\k)  can  itself  be  a  vector.  As  apparent,  vr(k)  is  constituted  by  the  samples  of 
the  original  signal  over  the  interval  [r,  r  +  T  —  1] ,  T  being  any  positive  integer.  For  a  proper 
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comprehension  of  the  way  in  which  signal  lifting  reflects  into  the  input /output  representation 
of  a  system,  we  start  with  by  considering  the  effect  of  sampling  in  the  ^-domain.  Then,  we 
will  pass  to  analyzing  the  action  of  lifting  on  time-invariant  and  T-periodic  systems,  both  in 
state-space  and  input/output  forms. 

5.2.1  Sampling  in  the  z-domain 

Again,  consider  a  discrete-time  signal  v(-),  whose  z- transform  is  given  by  the  celebrated 
formula  (written  in  terms  of  a  instead  of  z,  and  t  instead  of  k,  as  more  typically  done) 

OO 

V(a)  =  ^^v(£)cr_t  (5.44) 

t= o 

Hence,  the  ^-transform  of  the  periodically  sampled  signal 

v®(k)  =  v(kT  +  i-  1)  z  =  1,2,...  ,  T  (5.45) 

is  given  by 

OO 

V®(z)  =  ^Tv{kT  +  i~  l)z~k  (5.46) 

k= 0 

As  can  be  seen  in  the  expressions  above,  we  use  different  symbols  for  the  complex  variables 
entering  the  ^-transforms:  a  for  time  t  and  z  for  the  sampled  time  k.  In  fact,  we  are  using 
two  different  time  scales,  shown  in  Fig.  5.1.  In  the  discrete  time  scale  t  there  are  T  samples 
for  each  sample  of  the  time  scale  k.  In  other  words,  k  indicates  the  number  of  periods,  and 
t  the  number  of  samples,  with  T  samples  per  period,  i.e. ,  t  =  Tk.  The  operator  a  advances 
t  by  one,  while  the  operator  z  advances  k  by  one.  Because  it  takes  T  shifts  a  of  time  t  to 
make  one  shift  z  of  period  k,  we  have  z  =  aT. 

Now  denote  by  8(t)  the  Kronecker  symbol,  i.e.,  5(0)  =  1,  and  5{t)  =  0  for  t  ^  0.  With 
this  symbol  we  can  write 

OO 

v(9)  =^v(t)S(t-g)  (5.47) 

t= o 
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Figure  5.1:  Time  scales  used  in  lifted-time  reformulation;  period  T  is  equal  to  5. 


because  the  only  nonzero  term  in  the  summation  will  be  that  for  t  =  q,  because  only  for 
that  value  of  t  will  S(t  —  q)  be  equal  to  1.  Similarly 

OO 

v(q)a~q  =  ^  -  q)  (5.48) 

t= o 

and  if  we  let  q  =  kT  +  i  —  1 

OO 

v(kT  +  i  —  l)a~<ykT+l^l">  =  —  ( kT  +  i  —  1)]  (5.49) 

t= o 

or,  multiplying  both  sides  by  al~l 

OO 

v(kT  +  i  -  l)a~kT  =  -  kT  -  i  +  l)^"1  (5.50) 

t= o 
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Evaluating  Eq.  (5.46)  in  z  =  a7  it  can  be  seen  that 


where 


V^(aT 


a*  = 


^  v(kT  +  i—  l)c 


using  Eq.  (5.50) 

OO  CX) 

=  u(£)q~*£(t  —  kT  —  i  +  l)^1 


fc=0  4=0 


^^^^u(t)(T  t5(t  —  kT  —  i  +  1)  a*  1 


k= 0  4=0 


rearrange  some  terms  inside  the  square  brackets 

[  1 


OO  OO 


-  kT  -  i  +  1) 


_—t  i—  1 

a  a 


4=0  k= o 


=f  fi(t) 


^2v(t)fi(t)a  f  a1  1 


(5.51) 


fi(t)  =  '52s(t-kT~i+ 1) 


(5.52) 


is  a  T  periodic  sequence  of  discrete  pulses.  Inside  the  parenthesis, 

t  is  the  independent  variable  (will  change  with  outer  summation  index) 
kT  changes  with  the  summation  index 
i  changes  with  the  specific  choice  of  /)(f) 

Now,  recall  definition  (5.39)  of  the  T-roots  of  the  unit 


/ 2ni\  /2tt\  .  .  / 2tt 

=  exp  —  =  cos  —  +  »sm  — 


(5.39)  repeated 


and  the  associated  property  (5.40) 


1  s  =  0,  ±T,  ±2 T, . . . 
0  otherwise 


(5.40)  repeated 
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which  we  can  rewrite,  with  s  —  i  —  1  —  t 


1  y^1  Mi-i-t)  _  f  1  i-l-t  =  0,±T,±2T,... 

T  ^  1  0  otherwise 

fc=0  k 


Consider  the  conditions  under  which  the  summation  is  equal  to  zero: 


( 


< 


0 

±T 

_l_2 t  which  imply 


i  —  1  —  t  = 

i  —  1  —  t±T  = 
i  -  1  -  t  ±  2T  = 


and,  multiplying  by  -1 


t  —  %  +  1  =  0 

t  tT  -  i  +  1  =  0 

l  t  T  2T  —  i  +  1  =  0 


0 

0 

0 


(5.53) 


(5.54) 


(5.55) 


Comparing  these  conditions  with  those  under  which  the  summation  inside  Eq.  (5.52)  is  equal 
to  zero,  it  can  be  observed  that  they  are  the  same,  and  therefore  we  can  write 

1  T_1 

=  (5.56) 

k= 0 


Note  in  passing  that  this  expression  is  just  the  Fourier  sum  of  the  periodic  function 
Now  rewrite  Eqs.  (5.51)  and  (5.52): 


V^\aT) 


OO 

(*)/<(*)  O'-* 

_t= 0 
oo 


<7*_1 


(5.51)  repeated 


hit)  =  J2^t~kT-i  + 1) 


(5.52)  repeated 


k= 0 
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Use  Eq.  (5.56)  for  the  underlined  portion 


a  = 


VXt)Mt)a  * 


t= o 

oo 


a 


i- 1 


T—l 


-‘"‘I  •  a~' 


1= 0 


k= 0 


rearrange: 


rearrange  again: 


T 

1 

T 

1 

T 


T- 1 


cr 


i— 1 


1=0 


fc=0 


i-1 


T—l  oo 

:<-i) 


/c=0  1=0 

T—l 

^U(a0fe)0fc(i-1} 


— t 


a 


_i— 1 


a 


,k=0 


i-1 


To  understand  the  final  step  leading  to  Eq.  (5.57),  rewrite  Eq.  (5.44): 

OO 

rM  =  £  v(t)a  *  (5.44)  repeated 


1=0 


and  make  the  substitution  a  <—>  acj)k: 


(5.57) 


OO 

VW1)  =  £  v(t)(<r^)-‘ 

1=0 

Conversely,  the  2-transform  U(cr)  of  the  original  signal  u(t)  can  be  recovered  from  the 
2-transform  V^i)(z)  of  the  sampled  signals  v(kT  +  i).  To  see  this,  start  from  Eq.  (5.44): 

OO 

v(o)  =  £  v(t)a  *  (5.44)  repeated 

1=0 

rewrite  t  as  t  =  iT  +  k  (where  i  counts  the  number  of  periods  and  k  counts  the  sample  points 
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within  one  period) 


oo  T—  1 


a 


rearrange: 


add  one  to  k: 


EE  v(iT  +  k) 

i= 0  k= 0 

T- 1  oo 

J2<j-k^2v(iT  +  k ) 

A:=0  i=0 


-(iT+fe) 


(X 


-(iT) 


J2a~k+1J2v(iT + fc+1) 


a 


-(iT) 


(5.58) 


fc=l  i=0 


For  the  next  step  in  the  derivation  of  V(a),  first  rewrite  Eq.  (5.46): 


V^iz) 

OO 

=  J2v(kT  +  i-  l)z~k 

(5.46)  repeated 

change  notation  k  i: 

V^iz) 

k= 0 

oo 

=  J2v{iT  +  k-  1)^ 

i=0 

change  notation  again  z  aT : 

V^k\aT) 

oo 

=  vjiT  +  k  —  1  )cr~lT 

i= 0 

which  can  be  substituted  into  Eq.  (5.58)  to  get  the  converse  relation 

T 

V(a)  =  V{k\aT)a~k+l  (5.59) 

k=l 
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Chapter  6 


Constant-Coefficient  Representations  of  Periodic-Coefficient  Discrete  Linear 
Systems 

6.1  Introduction 

Various  types  of  active  control  of  a  helicopter  rotor  have  been  proposed  as  a  means  to  reduce 
vibration  and  noise,  and  to  improve  performance.  Very  popular  schemes  are  known  as  Higher 
Harmonic  Control  (HHC)  and  Individual  Blade  Control  (IBC),  typically  defined  as  active 
controls  implemented  in  the  fixed  system  and  in  the  rotating  system,  respectively.  Extensive 
theoretical  and  experimental  research  has  been  performed  on  HHC  and  IBC  for  nearly  three 
decades,  and  the  first  implementations  on  production  aircraft  are  currently  being  considered. 

Typical  HHC/IBC  implementations  are  based  on  frequency  domain  update  formulas  of 
the  general  type 

Zn+ 1  =  Zn  +  TnA6n  (6-1) 

where  Z  and  A 6  contain  harmonics  of  outputs  (usually  a  vector  of  hub  load  components)  and 
inputs  (a  vector  of  HHC/IBC  inputs),  T  is  a  transfer  (or  gain)  matrix  of  suitable  dimensions, 
which  depends  on  the  frequency  response  of  the  helicopter,  and  n  is  an  index  denoting  the 
rotor  revolution.  Both  the  stability  of  the  update  formula  of  Eq.  (6.1),  and  the  effects  on 
rotor  vibrations,  required  power,  and  noise,  have  been  studied  extensively. 

Although  in  steady  state  conditions  they  only  apply  harmonic  input  signals  to  the  rotor, 
HHC  and  IBC  are  closed-loop  dynamic  control  systems  that  act  by  modifying  the  rotor  dy¬ 
namics,  so  it  is  very  important  that  they  do  not  affect  negatively  the  aeromechanical  stability 
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of  the  helicopter.  As  for  every  active  rotor  control  system,  the  potential  for  aeroservoelas- 
tic  instabilities  exists,  especially  for  hingeless  and  bearingless  rotor  configurations,  which 
typically  have  lowly  damped  rotor  inplane  modes  or  coupled  rotor-fuselage  modes. 

However,  until  very  recently,  essentially  no  methodology  existed  in  the  published  litera¬ 
ture  to  study  the  coupled  rotor-fuselage  stability  of  a  helicopter  with  an  HHC  or  1BC  system 
turned  on,  and  no  quantitative  results  were  available.  One  probable  contributing  factor  is 
that  Eq.  (6.1)  does  not  readily  lend  itself  to  a  state-space  continuous  representation  that, 
when  coupled  to  linearized  models  of  rotor-fuselage  dynamics,  could  provide  closed-loop  sta¬ 
bility  information  in  the  form  of  stability  eigenvalues  or  Floquet  characteristic  exponents.  In 
fact,  a  state-space  representation  requires  that  it  be  possible  to  describe  the  entire  dynamics 
of  the  system  from  a  snapshot  at  any  instant  in  time,  and  there  are  at  least  two  aspects  of 
Eq.  (6.1)  that  appear  to  violate  this  requirement.  Firstly,  the  extraction  of  the  harmonics 
in  Z  requires  a  harmonic  analysis,  which  cannot  be  carried  out  at  a  single  instant  in  time, 
but  requires  a  time  interval,  of  at  least  one  cycle.  Secondly,  Eq.  (6.1)  describes  an  update 
at  discrete  intervals  in  time. 

The  situation  has  been  partially  remedied  by  Refs.  [9]  and  [28].  In  both  studies,  the  basic 
approach  was  to  convert  the  discrete  components  of  the  model  to  continuous,  or  to  neglect 
their  discrete  nature  altogether.  While  this  approach  represents  a  very  useful  first  step,  the 
system  described  by  Eq.  (6.1)  is  discrete  and  should  rigorously  be  treated  as  such. 

A  key  ingredient  of  a  discrete  stability  and  response  analysis  is  the  reformulation  of  a 
discrete  linear  system  with  periodic  coefficients  into  an  equivalent  one  with  constant  coeffi¬ 
cients.  Therefore,  the  main  objective  of  this  chapter  is  to  present  four  techniques  to  achieve 
this  goal,  namely:  (i)  time-lifted,  (ii)  cyclic,  (iii)  frequency-lifted,  and  (iv)  Fourier  reformu¬ 
lations.  The  theoretical  development  of  these  techniques  is  discussed  in  detail  in  Ref.  [4], 
This  chapter  summarizes  the  results  of  Ref.  [4]  of  potential  interest  for  the  helicopter  dy¬ 
namics  community,  and  presents  the  application  of  technique  (i)  to  a  simple  example,  i.e., 
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the  one-DOF  flapping  equation  for  an  isolated,  rigid,  articulated  rotor  blade  in  a  rotating 
coordinate  system. 

6.2  Discrete  formulation  of  the  blade  flapping  equation 


The  example  that  will  be  used  in  this  chapter  is  the  rigid  flapping  equation  of  motion  of  an 
isolated  rotor  blade  hinged  at  the  axis  of  rotation  (see,  for  example,  Ref.  [24]): 


fi  +  ^fl  +  ^/xsin  ip\  f3  +  1 +  ^  cos  ip  + /J?  sm2ip\  [3  = 

=  o  fl  +  yRsin-^  +  2/j2  sin2  ip)  9 (if})  -  ^  \  +  2/xsin'0^)  A  (6.2) 


which  can  be  rewritten  in  first-order  form  as: 

|  x  =  Ac(ip)x  +  Bc(ip)u 

\  y  =  Cc(ip)x  +  Dc(ip)u 

where  x  and  u  are,  respectively,  the  state  and  control  vector,  defined  as 


and  the  output  vector  y  is  chosen  to  be  equal  to  x.  The  system  matrices  are: 


(6.3) 


(6.4) 


A  .c(ip) 

BcW 


0  1 

A2l(lp)  -^22  (VO 

0  0  0 

f(ip)  f  ('ll))  cos  ip  f  (ip)  sin  ip 


(6.5) 

(6.6) 


with 


CM 

Dc(ip) 


0  0  0' 
0  0  0 


(6.7) 

(6.8) 


fW  =  l 


1  +  -/isin-0  +  2/j2  sin2  ip 
o 
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and 


^21  WO 
^22  WO 


1  + 


7/4 


-  /a  cos  ip  +  /i2  sin  2-0 


-  (  1  +  -/isin0 


Because  of  the  simplicity  of  the  flapping  equation,  all  the  system  matrices  above  could  be 
easily  obtained  by  inspection.  For  more  sophisticated  mathematical  models  the  matrices 
would  typically  be  obtained  numerically,  by  perturbing  the  governing  equations  about  a 
steady-state  equilibrium  position. 

In  a  discrete  representation,  the  system  is  defined  at  discrete  time  intervals  ipk  =  kApk  = 
1,2,...  .  The  corresponding  state-space  description  is  given  by 

(  x(/c  +  1)  =  A(k)'x.(k)  + 


(6.9) 


y(*0 


=  C(k)~x(k)  +  D(k)u(k) 


where,  for  example,  x(/c)  is  a  short-hand  notation  for  x(fcA0).  In  a  periodic  system,  such 
as  a  helicopter  rotor  blade,  all  the  matrices  have  a  common  period  equal  to  T  samples,  that 
is,  A{k)  =  A{k  +  £T)  with  £  =  1,  2,  3, . . .  Equation  (6.9)  yields  the  solution  x(/c  +  1)  at  time 
0  =  (k  +  1)A0  in  terms  of  the  state  x(/c)  and  controls  u (k)  at  time  0  =  kA'i/j. 

The  starting  point  to  link  the  continuous  system  matrices  of  Eq.  (6.3)  to  the  discrete 
system  matrices  of  Eq.  (6.9)  is  the  continuous  time  solution  formula  [34] 

x(l)  =  $01,r)x(r)  +  J  §A{t,(j)Bc(a)u(a)  da  (6.10) 


which  gives  the  solution  vector  x(t)  starting  from  the  initial  time  t  —  r,  with  an  input 
vector  u(f)  and  initial  conditions  x(r).  The  matrix  r)  is  the  state  transition  matrix 

of  the  system.  Therefore,  the  solution  x(/c  +  1)  at  time  ip  —  (k  +  1)A ip,  starting  from  time 
ip  =  kAip,  and  assuming  that  the  input  vector  u  is  held  constant  over  that  time  interval  at 
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a  value  u  =  u (k),  is  given  by: 


x[(k  +  l)A-0]  =  [(k  +  l)Apk  Atp]  x(kAip) 

p(h-\-  l)A-0 

+  /  $a[(^  +  l)A]/cr  ].Bc(cr )  da  ■  u(cr)  (6.11) 

J  kA'tp 

Comparing  Eq.  (6.11)  and  the  first  row  of  Eq.  (6.9)  it  can  be  seen  that: 


m 

B{k) 


<1^4  [{k  +  1)A ipk  Aip] 


p(k-\-l)Axp 

J  kA'ijj 


$A[(k  +  l)Aipj  ]Bc(a)  da 


y  (rad) 


(6.12) 

(6.13) 


Figure  6.1:  Continuous  and  discrete  solutions  of  flapping  equation. 


To  compute  A(k),  recall  that  for  the  special  case  r  =  0,  the  i-th  column  of  $a(T  0) 
can  be  obtained  by  integrating  from  r  to  f  the  homogeneous  system  x  =  Ac(ip)x.  with  an 
initial  condition  vector  x(0)  equal  to  zero  except  for  its  i-th  component  which  is  set  to  1. 
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Therefore,  first  compute  the  two  transition  matrices  ^A[kAip  0]  and  $^[(/c+l)A'0  0].  Then, 
the  property  of  transition  matrices  [34] 

<J>(t2,0)  =  <J>(t2,ti)$(ti,0)  (6.14) 

can  be  used  to  obtain 

[(fc  +  1)A#  AiP\  =  [{k  +  1)A^  0]  [kA^  0] 

=  A(k)  (6.15) 

The  discrete  control  matrix  B(k )  can  be  obtained  from  Eq.  (6.13)  numerically  using  a 

trapezoidal  rule,  which  gives 

m 

B(k)  S3  A<7j  Qj$A[(fc  +  l)Ajjj  i\Bc(<Ji)  (6.16) 

i=0 

where  m  is  the  number  of  intervals  in  which  the  sample  Aip  has  been  divided  for  the  purpose 
of  computing  the  integral,  A  a  =  A -fm  ,  o%  =  kAijj  +  iAa,  and  a  —  1/2  if  i  —  1  or  i  —  m 
and  a  =  1  otherwise.  The  transition  matrix  $^[(fc  +  VjA'iJxr  ,]  is  computed  as  described  in 
the  calculation  of  A(k ),  and  is  therefore  given  by: 

^[(H1)A#  i]  =  $A[(k  +  l)Ail>  0]$^h,0]  (6.17) 

The  continuous  and  discrete  solutions  of  the  flapping  equation  for  a  typical  case  are 
shown  in  Fig.  6.1.  The  numerical  details  are  provided  in  the  Appendix. 

6.3  Transfer  Functions 

The  output  x(/c)  at  a  certain  time  kA'ip  of  a  linear  time-periodic  system  can  be  written  as  a 
linear  combination  of  the  past  values  of  the  input  up  to  time  kAi/)\ 

x(/c)  =  M0(k)u(k)  +  Mi(k)u(k  —  1)  +  M2(k)u(k  —  2)  +  . . . 

OO 

=  ^2  Mi(k)n{t  -  k)  (6.18) 

3=0 
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The  coefficients  Mfk)  are  periodic  matrices  with  a  period  equal  to  T  samples,  known  as 
periodic  Markov  coefficients.  Physically,  the  j-th  column  of  the  Markov  coefficient  Mk_fk) 
represents  the  output  response  x(/c)  of  the  system  to  an  impulse  applied  at  time  kAf  to  the 
j-th  component  of  the  input  vector. 

Express  now  the  generic  time  point  k  in  terms  of  the  total  number  £  of  elapsed  periods 
(or  rotor  revolutions),  i.e. 

k  =  £T  +  s  (6.19) 

For  example,  if  k  =  24  and  a  rotor  revolution  is  discretized  into  T  =  10  samples,  then  £  =  2 
(two  complete  revolutions  have  taken  place)  and  s  =  4  (this  is  the  fourth  sample  in  the  third 
revolution).  With  this  in  mind,  it  can  be  shown  [4]  that  the  response  of  a  time-periodic 
system  at  a  generic  time  instant  k  can  be  written  as  the  Enite  sum  of  the  response  of  T 
time- invariant  systems,  indexed  in  the  integer  s.  Therefore,  it  is: 

T—l 

x(£T  +  s)  =  ^xl,s(£)  (6.20) 

i= 0 

where  xitS(£)  is  the  output  of  a  time- invariant  system  having  Mfs),  Mi+T(s),  Mi+2r(s), . . . 
as  Markov  parameters.  Each  of  these  time-invariant  systems  is  forced  only  once  per  period. 
An  important  reason  why  this  can  be  done  is  the  periodicity  of  the  Markov  parameters,  so 
that  Mi+T(s  +  T)  =  Mi+T(s),  Mi+2T(s  +  2 T)  =  Mi+2T(s ), - 

Equation  (6.20)  expresses  the  fundamental  concept  that  the  response  of  a  system  with 
periodic  coefficients  can  be  expressed  in  terms  of  the  responses  of  multiple  systems  with 
constant  coefficients.  This  is  the  basis  for  invariant,  or  constant-coefficient,  representations 
of  linear  periodic  systems. 

Introduce  now,  as  usual  in  the  representation  of  discrete  systems,  the  “one-shift  ahead 
in  time”  operator  z.  The  shift  is  of  T  samples  in  time,  not  just  one  sample,  because  we  are 
interested  in  describing  the  behavior  of  each  of  the  constant-coefficient  systems  in  which  we 
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have  decomposed  the  original  periodic-coefficient  system.  In  other  words,  z  shifts  from  a 
generic  sample  IT  +  s  to  (£  +  1  )T  +  s,  rather  than  £T  +  s  +  1. 

We  can  then  define  a  sampled  transfer  function  as: 

OO 

Hi(z,s)  =f  J2MeT+i(s)z~£  (6.21) 

t=o 

and  rewrite  the  system  response  as: 

x(£T  +  s)  =  H0(z,s)u(£T  +  s)  +  H1(z,s)u(£T+s-l)  + ... 

+Ht-2(z,  s)u(£T  +  s  —  T  +  2)  +  s)u(ydT  +  s  —  T  +  1)  (6.22) 

=  H0(z,  s)u.(£T  +  s)  +  Ht~i(z,  s)z~1u(y£T  +  s  +  1)  +  . . . 

+H2(z,  s)z~1\i(IT  +  s  +  T  —  2)  +  Hi(z,  s)z~1u(£T  +  s  +  T  —  1)  (6.23) 

Note  that  now  the  output  at  a  generic  time  sample  £T  +  s  is  expressed  entirely  in  terms 
of  inputs  u  over  one  period  T,  written  backward,  Eq.  (6.22),  or  forward,  Eq.  (6.23),  in  time. 

The  Markov  coefficients  in  Eqs.  (6.18)  and  (6.21)  can  be  expressed  in  terms  of  the  system 
matrices: 


M0{k) 

=  D(k ) 

(6.24) 

Mi(k) 

=  C(k)B(k-  1) 

(6.25) 

M2(k ) 

=  C(k)A(k~l)B(k-2) 

(6.26) 

M3(k) 

=  C(k)A(k-2)B(k-3) 

(6.27) 

or,  equivalently: 

=  D(k) 

=  C{k)¥A{k)<$>A{k,k-i  +  l)B(k-i) 


M0(k ) 
M(.T+i 
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(6.28) 

(6.29) 


Substituting  into  Eq.  (6.21)  yields  for  the  transfer  functions  H(z,k),  after  some  manipula¬ 
tions: 


H0(z,k)  =  D(k)  +  C(k)[zI  -  ^(A;)]"1  $A(k,k-T  +  l)B(k)  (6.30) 

Hi(z,k)  =  zC(k)[zI  -^A(k)}-l<$>A(k,k-i  +  l)B(k-i) 

i  =  1,2,...  ,T-  1  (6.31) 

Additional  details  and  numerical  values  for  the  flapping  blade  example  can  be  found  in  the 
Appendix. 

A  second  type  of  transfer  function  is  very  useful  for  the  construction  of  time-invariant 
representations  of  periodic  systems.  Using  the  unit  delay  operator  a-1  (e.g.,  u(k  —  1)  = 
cr_1u(A;)),  Eq.  (6.18)  can  be  rewritten  as 

x(/c)  =  [M0(k)  +  +  M2{k)a~2  +  . . .  ]  u (k) 

=  G(a, -)\k  u(k)  (6.32) 

where  the  input- output  transfer  operator  G (a,  • )  | is  defined  as 

G(a,  -)|fc  =  [. M0(k )  +  M\{k)a~l  +  M2(k)a~2  +  . . .  ]  (6.33) 

Note  that  G(a,  -)|fc  in  Eq.  (6.32)  is  an  operator  on  u (k),  not  a  multiplicative  factor  for  u (k), 
i.e. ,  it  is  a  more  elaborate  way  of  writing  f(u(k)). 

If  the  input-output  transfer  operator  is  evaluated  in  a  specific  time  point,  say  k,  it  results 
in  a  transfer  function,  denoted  by  G(a,k).  Because  they  both  describe  the  same  dynamic 
system,  this  transfer  function  is  obviously  related  to  the  sampled  transfer  function  Hi(z,s), 
defined  in  Eq.  (6.21).  It  can  be  shown  [4]  that  the  two  transfer  functions  are  related  by: 

T— 1 

G{a,  k)  =  Y l  (6.34) 

i= 0 
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and  the  inverse  relationship 


i  r-1  .1 

Hi(aT,k)  =  —  k)4>h  a1  (6.35) 

U=o 

where  (f)  =  exp(27r j/T)  =  cos(27 r/T)  +  j  sin(27r/T).  Note  that  if  the  original  periodic  system 
is  time-invariant,  i.e.,  T  =  1  and  hence  a  —  z,  then  the  transfer  function  G(a,k )  is  not 
a  function  of  k  and  reduces  to  the  transfer  function  G(a)  of  the  time-invariant  system. 
Additional  details  and  numerical  values  for  the  flapping  blade  example  can  be  found  in  the 
Appendix. 

6.4  Time-invariant  reformulations 
6.4.1  Time-lifted  reformulation 

The  operation  of  “lifting”  consists  in  “packing”  the  values  of  a  signal  over  one  period  into  a 
new  enlarged  signal.  Time-lifting  results  in  a  constant  coefficient  state-space  system  of  the 
type 

x(fc)(£  +  l)  =  Fk^k\£)  +  Gkuk(£)  (6.36) 

y(£)  =  Hk^k\£)  +  Ekuk(£)  (6.37) 

where  x^(£)  is  a  shorthand  notation  for  x(£T  +  k),  with  £  the  number  of  periods  elapsed 

from  time  zero,  and  k  the  time  within  the  period  at  which  the  sample  is  taken.  The  vector 

uk(£)  is  the  “packed”  input  vector,  defined  as: 

u(£T  +  k) 
u(£T  +  k  +  l) 
u(£T  +  k  +  2) 

u {IT  +  k  +  T  -  2) 
u  (£T  +  k  +  T  -  1) 

The  packed  output  vector  yk(£)  is  defined  in  a  similar  way.  Therefore,  if  the  system  has  n 
states,  m  controls,  and  p  outputs,  and  is  sampled  T  times  over  one  period,  the  dimensions 


(6.38) 
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of  the  state  matrix  Fk,  the  control  matrix  Gk,  the  output  matrix  Hk ,  and  the  feed-through 
matrix  Ek,  will  be  n  by  n,  n  by  mT,  pT  by  n,  and  pT  by  mT,  respectively. 

The  matrices  in  Eqs.  (6.36)  and  (6.37)  are  given  by  [4]: 


Fk  =  (k) 


Gk  = 


with: 


G 


k  1 


G 


k2 


G 


kT 


(6.39) 

(6.40) 


Hk  = 


Gkl  =  $A(k  +  T,k  +  l)B(k ) 

Gk2  =  ®a(^  +  T,  k  +  2)B(k  +  1) 

GkT  =  $A(k  +  T,k  +  T)B(k  +  T-  1) 
C(k) 

C(k  +  l)$A(k+l,k) 


Ek 


C{k  +  T  -  l)<S>A{k  +  T  -  1,  k) 
0 

D[k  +  i  —  1) 


(6.41) 


i  <  ] 

{Ek)a  =  {  C(k  +  i-l)<f>A(k  +  i-l,k  +  j)x 

xB(k  +  j  —  1)  i  >  j 


(6.42) 


For  the  derivations  of  the  present  chapter  it  will  be  £  =  0,  i.e.,  the  starting  point  is  the  first 
point  of  the  period.  Also,  the  subscript  k  will  be  dropped  from  all  the  system  matrices. 

The  time-lifted  matrix  F  is  given  by  Eq.  (6.39)  as  F  —  Ta(0),  where  Ta(0)  '=  $a(T,  0) 
is  the  transition  matrix  at  the  end  of  one  period.  Therefore,  the  periodic  system  will  be 
stable  if  and  only  if  the  time-lifted,  constant  coefficient  system  is  stable.  The  time-lifted 
state  matrix  F  is  given  by  ^[lOA-^  0]  =  $ A [6. 2832, 0],  in  the  Appendix. 

Note  that  the  open  loop  stability  characteristics  of  the  discrete  systems,  in  both  the  peri¬ 
odic  and  the  time-lifted  formulations,  are  identical  to  those  of  the  corresponding  continuous 
system.  In  fact,  the  discrete  state  matrices  A(k),  and  therefore  also  dy^O),  have  been  built 
from  the  continuous  time  solution  formula,  Eq.  (6.10).  Therefore  the  transition  matrices  at 
the  end  of  one  period  for  the  continuous  and  the  discrete  systems  are  identical  and  identical 
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are  their  eigenvalues,  which  determine  the  stability  levels. 


From  Eq.  (6.40),  the  time-lifted  control  matrix  G  for  the  flapping  blade  example  is  given 


by: 


G  = 


G\ 


G 5 


with: 


G 


T—l 


(6.43) 


G1  =  <S>A(T,1)B(0) 

G2  =  $a(T,2)B(1) 

GT-i  —  $a(T,  T  —  1)B(T  —  2) 

Gt  =  ^a(T,T)B(T-  1) 

The  control  matrices  B(. . .)  are  given  in  Table  in  the  Appendix.  The  transition  matrices 
&a(T,  1),  <&a{T,2),  . . .  can  be  obtained  from  the  transition  matrix  property,  Eq.  (6.14)  as 
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follows 


*a{T,T) 


®a{T,T-  1) 


$a{T,T-  2) 


*a{T,T-  3) 


A(T  -  1) 

<f>A(T,T-l)$A(T-l,T-2) 
<Fa(T,T-  l)A(T-2) 

A(T -*  X)A(T  —  2) 

<f>A(T,  T  —  2)$a(T  —  2,  T  —  3) 
$a(T,  T  —  2)A(T  —  3) 

,4(7'  -  l),4(7'-2),4(7'-3) 


<f>A(T,l)  =  $a(T,  2)$a(2,  1) 

=  $a(T,2)A(1) 

=  A(T  —  1)A(T  —  2)A(T  —  3)  x 
...  x  A(2)A(1) 


where  the  state  matrices  A(. . . )  are  given  in  the  Appendix.  For  the  blade  flapping  equation, 
which  has  2  states  and  3  controls,  and  for  a  period  discretized  in  T  —  10  intervals,  the 
time-lifted  control  matrix  G  has  2  rows  and  30  columns  (see  the  Appendix,  Eq.  (6.70)  for 
the  numeric  values  of  the  matrix). 

From  Eq.  (6.41),  the  time-lifted  output  matrix  H  for  the  flapping  blade  example  is  given 


by: 


C(  0) 

I 

C(1)$A(1,0) 

<Mi,o) 

H  = 

C(2)$a(2,0) 

= 

$a(2,0) 

C(9)$a(9,0) 

<M9,0) 
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because  for  the  blade  example  all  the  discrete  output  matrices  C(k)  are  taken  to  be  equal 
to  the  identity  matrix.  With  2  states  and  2  outputs,  and  for  a  period  discretized  in  T  =  10 
intervals,  the  time-lifted  output  matrix  G  has  20  rows  and  2  columns.  The  transition  matrices 
&A(k,  0)  and  the  output  matrix  H  are  given  in  the  Appendix. 

Finally,  from  Eq.  (6.41),  the  time-lifted  feedthrough  matrix  E  for  the  flapping  blade 
example  is  given  by: 


E 


Eij  — 


02x3  i  <  j 

D{i  -  1)2x3  =  0  i  =  j 

=  $A(i  -  hj)B(j  ~  1)  i  >  3 


(6.45) 


because  for  the  blade  example  all  the  discrete  output  matrices  C(k)  are  taken  to  be  equal  to 
the  identity  matrix  and  all  the  D(k)  are  zero  matrices.  The  transition  matrices  &A(i  —  1  ,j) 
are  obtained  using  <&A(i  —  1  ,j)  =  <&A{i  —  1 , 0) (j,  0) .  The  notation  E^  indicates  the  i-th 
2x3  block  rowwise  and  j-th  columnwise  of  E.  With  2  outputs  and  3  controls,  and  for  a 
period  discretized  in  T  =  10  intervals,  the  time-lifted  feedthrough  matrix  E  has  20  rows  and 
30  columns.  The  matrix  is  given  by  Eq.  (6.72)  in  the  Appendix. 


6.4.2  Cycled  reformulation 


In  the  time-lifted  reformulation  shown  in  the  previous  section,  the  state  vector  remains  of  its 
original  dimension,  i.e.,  it  is  not  lifted  or  augmented  in  any  form.  In  the  cycled  reformulation, 
all  the  state,  input,  and  output  vectors  are  augmented.  For  a  generic  signal  vector  v(fc), 
which  can  represent  any  of  those  vectors,  the  operation  of  “cycling”  results  in  a  vector  Vr\k) 
defined  as  [4]: 


y(i) 


(k) 


v(k)  k  =  t  +  i  —  1  +  tT 
0  otherwise 


(6.46) 


where  the  superscript  ( % )  in  Vr\k)  denotes  the  partition  of  the  cycled  vector  vT(k).  In  other 
words,  at  each  time  point,  the  vector  vT(k)  has  a  unique  nonzero  sub- vector,  which  cyclically 
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shifts,  e.g. 


vT(r)  = 

1 

< 

... 

_ 1 

,  VT(r  +  1)  = 

0 

v(t  +  1) 

0 

0 

0 
0 

,•••  ,  vt(t  +  1)  = 

I 

v(t  +  T  - 

Cycling  results  in  a  constant  coefficient  state-space  system  of  the  type 


(6.47) 


xr(/c  +  1) 

y(k) 


FT-kT{k)  +  GTuT(k ) 
HTxT(k)  +  ETuT(k) 


(6.48) 

(6.49) 


where  the  “state”  matrix  Fr  is  given  by: 


Ft  = 


0 

0 

0 

A{t  +  T  —  1) 

A(t) 

0 

0 

0 

0 

A(t  +  1) 

0 

0 

(6.50) 

0 

0 

. ..  A(r  +  T-  2) 

0 

the  control  matrix  Gr  has  exactly  the  same  structure  as  FT,  but  with  B(-)  instead  of  A(-), 


the  output  matrix  Hr  is  given  by: 

'  C{t)  0 

0 

0  C(t  +  1)  . . 

0 

(6.51) 

Hr  = 

0  0 

.  C(t  +  T-  1) 

and  the  feedthrough  matrix  Er  has  exactly  the  same  structure  as  FtT1  but  with  D(-)  instead 

ofC(-). 

Considering  the  matrix  Fr,  it  can  be  shown  [4]  that  its  eigenvalues  coincide  with  all 
the  T-th  roots  of  the  characteristic  multipliers  of  the  discrete  state  matrix  A(-).  Therefore, 
the  periodic  system  is  stable  if  and  only  if  its  cycled  reformulation  is  stable.  Additional 
properties,  important  from  a  control  system  standpoint,  can  be  found  in  Ref.  [4], 
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6.4.3  Frequency-lifted  reformulation 


The  frequency  lifting  is  conceptually  similar  to  the  time  lifting,  except  that  in  this  case  the 
enlarged  vectors  contain  all  the  harmonics  of  the  input  and  of  the  output  signals.  Therefore, 
this  is  a  frequency  domain  representation.  This  reformulation  is  naturally  set  in  the  ^-domain 
directly. 

Consider  a  generic  discrete  signal  v(fc),  the  ^-transform  of  which  is  V(cr)  (recall  that 
while  z  shifts  signals  by  one  sample,  a  shifts  signals  by  T  samples,  and  therefore  it  is  the  z- 
transform  equivalent  when  we  wish  to  consider  the  behavior  of  the  system  every  T  samples). 
Then,  we  can  define  a  frequency-lifted  vector  V  (a)  as: 


V(a) 

VM) 

V(a02) 


> 


(6.52) 


(  V(a0T-1)  J 

where  again  <f>  =  exp(j'27r/T).  Notice  that  the  i-th  component  of  this  vector  is  nothing  but 
the  ^-transform  of  the  complex  signal  obtained  by  multiplying  the  signal  v(f)  by  the  complex 
number  cj)^~l+1lk. 


When  the  general  definition  of  Eq.  (6.52)  is  specialized  to  the  state  vector  x(fc),  the 
control  vector  u(k),  and  the  output  vector  y(k),  we  obtain  frequency-lifted  vectors  X(cr), 
U(cr),  and  Y(cr),  each  with  T  subvectors,  and  total  size  of  nT,  mT ,  and  pT,  respectively. 
The  frequency-lifted  input  vector  U(cr)  and  output  vector  Y(cr)  are  related  by: 


Y  (cr)  =  W  (cr)U(cr) 


(6.53) 


where  W (a)  is  the  frequency-lifted  transfer  function  of  the  periodic  system.  The  matrix  has 
dimensions  pT  by  mT,  and  it  is  composed  of  T  by  T  blocks  of  p  rows  and  m  columns.  The 


submatrix 


W(a) 


<j+l,r+l 


in  row  q  +  1  and  column  r  +  1,  with  r,  q  —  0, 1,  2, ...,  T  —  1  is  given 
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(6.54) 


by  [4] 


Therefore,  for  example 


W(a) 


<2+1, r+1 


1  T-1 

-ycK.v'1'-'1 

£=0 


q  =  0,  r  =  0 
q  =  1,  r  =  0 
q  =  0,  r  =  1 


W(a) 


T- 1 


T 

1 

T 

1 

T 


£=0 
T— 1 


^=o 

T—l 


£=0 


All  these  matrices  will  generally  have  complex  coefficients.  Note  that  the  various  elements 
of  the  frequency-lifted  transfer  function  are  expressed  in  terms  of  the  input-output  transfer 
function  G(a,  k). 


6.4.4  Fourier  reformulation 


The  Fourier  reformulation  is  closely  related  to  the  frequency-lifted  reformulation,  and  is  also 
a  frequency-domain  representation.  This  formulation  is  obtained  by  considering  the  periodic 
system  in  an  Exponentially  Modulated  Periodic  (EMP)  regime. 

A  signal  v(/c)  is  said  to  be  EMP  if  there  exists  a  (complex)  number  A  ^  0  such  that 

v(£T  +  s)  —  v(s)  XeT  (6.55) 


where  s  is  any  of  the  T  samples  that  make  up  one  period.  For  A  =  1  or  any  T-root  of  1  an 
EMP  signal  is  also  periodic.  Any  EMP  signal  v(fc)  can  also  be  written  as  the  following  EMP 
Fourier  expansion  [4] 

T—l 

v(k)  =  ^^V£(j)ekXk  (6.56) 

£=0 

where  \q  are  the  coefficients  of  the  Fourier  expansion  of  the  periodic  signal  v(fc)  =  v(fc) X~k. 
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To  obtain  the  Fourier  reformulation,  the  state  vector  x(fc),  the  control  vector  u (k),  and 
the  output  vector  y (k),  are  all  expressed  as  EMP  Fourier  expansions.  For  example 

T—l 

u (k)  =  uiaekak  (6.57) 

e=o 

(see  Ref.  [4]  for  details  of  the  substitution  of  Xk  in  Eq.  (6.56)  with  ak  in  Eq.  (6.57)).  The 
system  matrices  are  similarly  expanded.  For  example,  for  the  state  matrix  A(k)  we  have 

T—l 

A(k)  =  J2  M£k  (6.58) 

e=o 

Substituting  the  EMP  expansions  of  the  system  matrices  into  the  system  equations, 
Eq.  (6.9),  and  equating  all  terms  at  the  same  frequency,  we  obtain  the  following  matrix 
equation 


aj\fk  =  Ak  +  Bvl 


(6.59) 


y  =  Ck  +  Vu 


(6.60) 


where  x,  u,  and  y  are  vectors  formed  with  the  harmonics  of  x(k),  u (k),  and  y(k),  respectively, 
organized  in  the  following  fashion 


x  =  < 


Xo 

Xl 

XT-1 


(6.61) 


and  similarly  for  u  and  y.  A ,  B ,  C,  and  V  are  formed  with  the  harmonics  of  A(k),  B(k),  C(k), 
and  D{k ),  respectively.  For  example 


Aq  At-  i  At~2  •••  Ai 

Ai  A0  At-i  ...  A2 

A=  :  :  :  j 

At—2  At~  3  At—4  ■  ■  ■  AT-i 

At-i  At -2  At -3  ■  ■  ■  A0 


(6.62) 
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and  similarly  for  B,  C ,  and  T>.  The  matrix  J\f  is  the  block  diagonal  matrix 

AT  =  blkdiag  [<$>k  I]  k  =  0, 1,  2, T  -  1  (6.63) 

Then,  one  can  define  the  harmonic  transfer  function  as 

G{<j)  =  C[aNl  -  A]~lB  +  V  (6.64) 

This  transfer  function  provides  a  connection  between  the  input  harmonics  U  and  the  output 
harmonics  Y.  Note  that  for  the  steady-state  periodic  case  a  —  1.  Finally,  it  can  be  shown  [4] 
that  the  harmonic  transfer  function,  Eq.  (6.64),  coincides  with  the  frequency-lifted  transfer 
function,  Eq.  (6.53). 

6.5  Other  considerations 

Taken  together,  the  four  constant  coefficient  reformulations  described  in  the  previous  sections 
provide  the  foundation  required  to  study  a  coupled  rotor-fuselage  system  with  a  closed 
loop  1BC  or  HHC  system  entirely  as  a  discrete  system.  This  is  important  because  the 
effects  on  stability  and  performance  characteristics  of  discrete  elements  (e.g.,  control  updates, 
harmonic  analyses,  etc.)  and  sensor  and  computation  delays  can  then  be  studied  rigorously. 
The  approximations  intrinsic  in  converting  discrete  elements  to  continuous,  especially  if  the 
sampling  or  update  rates  are  low,  are  eliminated. 

Because  the  reformulated  systems  have  constant  coefficients,  all  the  control  system  de¬ 
sign  tools  for  discrete  time-invariant  systems  can  be  used  to  design  both  time-invariant  and 
periodic  rotor  controllers.  Possible  interactions  with  flight  control  systems  can  also  be  stud¬ 
ied,  and  gain  and  phase  margins  can  be  established.  However,  it  is  important  to  point  out 
that  while  a  periodic  coefficient  system  has  always  an  equivalent  constant  coefficient  rep¬ 
resentation,  the  reverse  is  not  true.  For  example,  a  constant  coefficient  system  will  have 
a  periodic  coefficient  analogue  only  if  the  direct  feed-through  matrix  Er  of  the  time-lifted 
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system,  Eq.  (6.45),  is  block  lower  triangular  [12].  This  has  practical  implications  on  the 
design  of  active  rotor  controls  in  the  rotating  system. 

6.6  Summary 

The  fully  discrete  modeling  of  a  coupled  rotor-fuselage  system  allows  a  rigorous  treatment 
of  closed-loop  HHC  or  1BC  systems,  including  the  discrete  elements  with  multiple  sampling 
and  update  rates. 

The  extraction  of  a  linearized  discrete  model  from  the  continuous  one  is  very  simple. 
The  Floquet  characteristic  multipliers  and  characteristic  exponents  are  the  same  for  the 
discrete  and  continuous  systems,  and  therefore  the  open  loop  stability  characteristics  of  the 
two  systems  are  exactly  the  same. 

Time-lifting  converts  a  system  with  periodic  coefficients  into  a  larger  system,  but  with 
constant  coefficients.  Time-lifting  is  based  on  two  properties  of  discrete  periodic  systems, 
namely:  (i)  the  response  of  the  periodic  system  can  be  modeled  as  the  combination  of  the 
responses  of  multiple  constant-coefficient  systems,  one  per  sample,  and  (ii)  using  shifting 
operators,  it  is  possible  to  “move”  all  the  system  inputs  at  the  same  sample  in  multiple 
periods  to  that  sample  in  a  single  period. 

Other  discrete,  constant  coefficient  reformulations  are  available,  and  can  be  useful  for 
theoretical  studies  and  practical  applications. 

6.7  Appendix 

This  appendix  contains  numerical  details  for  the  blade  flapping  examples  presented  in  this 
chapter.  The  parameters  used  in  the  flapping  equations  were:  Lock  number  7  =  8,  inflow 
ratio  A  =  0,  advance  ratio  /i  =  0.3,  collective  pitch  90  =  10  deg,  lateral  cyclic  6lc  =  2  deg, 
longitudinal  cyclic  6\s  =  3.  The  choice  of  setting  A  =  0,  not  very  realistic  from  a  practical 
point  of  view,  was  simply  made  to  simplify  the  calculations:  in  fact,  for  the  mathematical 
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model  of  this  study  A  is  neither  a  state  nor  a  control.  In  a  more  accurate  mathematical 
model,  A  would  probably  be  a  state,  have  its  governing  equations,  and  be  part  of  the  state 
vector  x.  To  take  it  into  account  in  the  present  study,  it  could  be  added  to  the  u  vector  as 
if  it  were  a  pseudo-control. 

The  transition  matrices  $a[^A ip  0],  k  =  1,2,...  ,10,  which  are  2  by  2  matrices,  are 
listed  below.  Note  that  $a[10A^  0]  is  the  transition  matrix  at  the  end  of  one  period,  the 
eigenvalues  of  which  determine  the  stability  of  the  periodic  system. 


$a[A^  0]  =  $^[0.6283,0] 


0.7849  0.4154 

-0.5769  0.3244 


$a[2A^  0]  =  $a[1.2566,0] 


0.4025  0.4616 

-0.5661  -0.1129 


$a[3A^  0]  =  $a[1.8850,0] 


0.1189  0.3466 

-0.3322  -0.2118 


$a[4A^  0]  =  $a[2.5133,0] 


-0.0285  0.2234 

-0.1538  -0.1730 


$a[5A^  0]  =  $a[3. 1416,0] 


-0.0919  0.1293 

-0.0578  -0.1299 


$a[6A^  0]  =  $a[3. 7699,0] 


-0.1086  0.0564 

0.0016  -0.1037 


$a[7A^  0]  =  $a[4.3982,0] 


-0.0929  -0.0011 
0.0460  -0.0779 


$a[8A^  0]  =  $a[5.0265,0] 


-0.0550  -0.0390 
0.0701  -0.0408 


$a[9A^  0]  =  $a[5.6549,0] 


-0.0113  -0.0512 
0.0637  0.0015 


$a[  10A^  0]  =  $a [6.2832,0] 


0.0195  -0.0399 
0.0319  0.0304 
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The  discrete  state  transition  matrices  A(k ),  k  —  1,2,...  ,10,  which  are  2  by  2  matrices, 


are  listed  below: 


7L(1)  =  7L(0.6283) 
A(  2)  =  4(1.2566) 
4(3)  =  A(1.8850) 
4(4)  =  4(2.5133) 
4(5)  =  4(3.1416) 
4(6)  =  4(3.7699) 
4(7)  =  4(4.3982) 
4(8)  =  4(5.0265) 
4(9)  =  4(5.6549) 
4(10)  =  4(6.2832) 


0.7849  0.4154  ' 

-0.5769  0.3244 

0.8030  0.3947  ' 

-0.5033  0.2965 

0.8467  0.3920  ' 

-0.3818  0.3153 

0.8923  0.4053  ' 

-0.2768  0.3638 

0.9108  0.4289  ' 

-0.2536  0.4232 

0.8946  0.4561  ' 

-0.3194  0.4809 

0.8622  0.4793  ' 

-0.4163  0.5244 

0.8346  0.4888  ' 

-0.4911  0.5309 

0.8134  0.4771  ' 

-0.5435  0.4823 

0.7936  0.4477  ' 

-0.5817  0.3974 


The  discrete  control  matrices  B(k),  k  =  1,2, . . .  ,10,  are  listed  below: 
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B(l)  = 

0.1792 

0.1720 

0.0419 

0.5319 

0.4847 

0.1996 

B(  2)  = 

'  0.2482 

0.1597 

0.1861  ' 

0.6993 

0.3643 

0.5848 

B(  3)  = 

'  0.2860 

0.0251 

0.2816  ' 

0.7739 

-0.0442 

0.7608 

B(  4)  = 

'  0.2700 

-0.1350 

0.2304  ' 

0.7108 

-0.4336 

0.5490 

B(  6)  = 

'  0.2085 

-0.1882 

0.0844  ' 

0.5394 

-0.5082 

0.1527 

B(  6)  = 

'  0.1360 

-0.1318 

-0.0268 

0.3517 

-0.3279 

-0.1103 

B(  7)  = 

'  0.0850 

-0.0567 

-0.0621 

0.2279 

-0.1294 

-0.1830 

B(  8)  = 

'  0.0654 

-0.0064 

-0.0643 

0.1889 

0.0062 

-0.1857 

B(  9)  = 

'  0.0745 

0.0382 

-0.0630  ' 

0.2298 

0.1433 

-0.1750 

B(10)  = 

'  0.1136 

0.1039 

-0.0426  ' 

0.3522 

0.3369 

-0.0829 

Transfer  functions  Hi(z,k )  and  G(a,k ) 

For  the  blade  flapping  example,  there  are  100  transfer  function  matrices  Hi(z,k)  because 
both  i  =  0, 1,  2, . . .  ,  T  —  1  and  k  =  1,  2, . . .  ,  T.  Each  of  the  matrices  has  2  rows  and  3 
columns.  Denoting,  for  convenience,  B'[k )  =  ^A(k,  k  —  T  +  l)B{k),  and  recalling  that 
D{k )  =  0,  the  expression  for  H0(z,  k ),  Eq.  (6.30)  becomes 


H0(z,k)=C(k)  [zI-HA{k)rlB\k) 


(6.65) 


which  represents  a  conversion  from  state  space  to  transfer  function  that  can  be  carried  out, 
for  example,  using  the  Matlab  function  ss2tf .  For  illustration,  consider  the  case  k  =  0,  i.e., 
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compute  Hq(z,0).  We  have: 


T4(0)  =  A(9)7L(8)  x  ...  x  A(0)  (6.66) 

$A(0,  -9)  =  A(9)A(8)  x  ...  x  7L(1)  (6.67) 


B( 0)  =  5(10),  (7(0)  is  the  identity  matrix,  and  5(0)  is  zero.  Using  ss2tf,  the  first 
column  of  Hq(z,  0)  can  be  computed  using  ss2tf  (A9*A8*A7*A6*A5*A4*A3*A2*A1*A0  and 
A9*A8*A7*A6*A5*A4*A3*A2*A1*B0,C,D, 1) ,  where  A0,  Al,  ...  denote  A(0),  A(l), ...,  re¬ 
spectively.  To  compute  the  second  and  third  columns,  the  last  argument  1  is  replaced  by  2 
and  3,  respectively. 

The  result  is 


H0(z,  0) 

x 


1 

U2  -  0.0499z  +  0.0019  X 
-0.0193^  +  0.0004  -0.0186^  +  0.0004 
0.0077^-0.011  0.0071* -0.011 


-0.0047;;  +  0.0001 
0.0028;;  -  0.003 


(6.68) 


Note  that,  although  the  100  Hi(z,k)  matrices  will  generally  have  different  coefficients, 
the  characteristic  polynomial  z2  —  0.0499;;  +  0.0019  will  be  the  same  for  all  of  them. 

The  transfer  function  G(a,k )  is  calculated  using  Eq.  (6.34): 

T-l 

G(cr,k )  =  Hi(aT ,  k)a~l  Eq.  (6.34)  repeated 
i= o 

For  the  blade  flapping  example  there  are  10  different  G(a,k ),  i.e.,  one  for  each  value  of 
k  —  0, 1, 2, . . .  ,  T  —  1.  For  example,  for  k  =  0  and  T  =  10 

9 

G(a,  0)  =  Hi(a9i  °)a"*  (6.69) 

i= 0 

The  first  term  in  the  summation,  i  —  0,  is  77o(<r9,0).  The  matrix  Hq(z,0)  is  given  by 
Eq.  (6.68).  Therefore  5o((J9,0)  is  obtained  from  Hq(z,0)  by  replacing  z  with  a9  in  all  the 
elements.  This  can  be  done  with  either  numerical  or  symbolic  manipulation  software. 
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Time-lifted  reformulation 


The  time-lifted  state  matrix  F  is  given  by  ^[lOA-^  0].  The  time-lifted  control  matrix  G  is 
given  by: 


G 


G\ 

G-2 

g3 

g4 

G5 

G6 

G7 

Gs 


Uri  w  2  ^3  W  ^5  ^6  ^7  ^8 


where: 

'  -0.0200  -0.0189 
-0.0064  -0.0065 

'  -0.0284  -0.0101 
-0.0077  -0.0052 

'  -0.0243  0.0138 

-0.0726  -0.2159 

0.1103  -0.0677 
-0.3155  0.1831 

"  -0.1904  0.0597 

0.2755  -0.0958 

'  -0.0675  0.1934 

0.0459  -0.0449 

'  0.1469  -0.0023 
0.0682  0.0064 

'  0.0382  -0.0630  ' 
0.1433  -0.1750 


0.0056 

-0.0306 

0.0004 

-0.0092 

-0.0312 

-0.0187 

-0.0869 

-0.0453 

-0.0038 

0.0136 

0.0064 

-0.2123 

0.0849 

0.2028 

-0.2506 

-0.2965 

0.2280 

-0.2139 

-0.1763 

0.1674 

-0.1142 

-0.1522 

0.0310 

0.0321 

0.1444  0.0745 
0.0670  0.2298 


(6.70) 
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The  time-lifted  output  matrix  H  is  given  by: 


Ht 


Hr 

H2 

h3 

h4 

h5 


[Hr  H2 

CO 

£ 

H 5  ] 

where: 

'  1.0000 

0  0.7849  -0.5769  ' 

0 

1.0000  0.4154  0.3244 

'  0.4025 

-0.5661 

0.1189  - 

0.3322  ' 

0.4616 

-0.1129 

0.3466  - 

0.2118 

"  -0.0285 

-0.1538 

-0.0919 

-0.0578 

0.2234 

-0.1730 

0.1293 

-0.1299 

'  -0.1086 

0.0016 

-0.0929 

0.0460 

0.0564 

-0.1037 

-0.0011 

-0.0779 

"  -0.0550 

0.0701 

-0.0113 

0.0637  ' 

-0.0390 

-0.0408 

-0.0512 

0.0015 

(6.71) 
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Finally  the  time-lifted  feedthrough  matrix  E  is  a  matrix  with  30  rows  and  20  columns. 
For  reasons  of  space,  only  the  first  4  rows  will  be  listed  below,  transposed: 


0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.1136 

0.1039 

-0.0426 

0.0000 

0.3522 

0.3369 

-0.0829 

0.0000 

0.2303 

0.2164 

-0.0669 

0.1792 

0.0472 

0.0476 

-0.0032 

0.5319 

0.2135 

0.2019 

-0.0579 

0.3602 

-0.0730 

-0.0676 

0.0246 

0.0993 

0.1609 

0.1527 

-0.0417 

0.3617 

-0.0857 

-0.0805 

0.0250 

-0.0636 

0.1098 

0.1046 

-0.0273 

0.3021 

-0.0771 

-0.0728 

0.0211 

-0.1186 

0.0631 

0.0604 

-0.0148 

0.2162 

-0.0721 

-0.0684 

0.0189 

-0.1535 

0.0198 

0.0193 

-0.0037 

0.1128 

-0.0641 

-0.0610 

0.0160 

-0.1705 

-0.0148 

-0.0138 

0.0048 

0.0108 

-0.0438 

-0.0418 

0.0103 

-0.1459 

-0.0329 

-0.0311 

0.0088 

-0.0608 

-0.0131 

-0.0127 

0.0024 

-0.0762 
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Chapter  7 


Discrete-Time,  Closed-Loop  Aeromechanical  Stability  Analysis  of  Helicopters 
With  Higher  Harmonic  Control 

7.1  Introduction 

Higher  Harmonic  Control  (HHC)  is  an  active  vibration  control  technology  that  tries  to 
reduce  N/rev  (N  being  the  number  of  rotor  blades)  vibratory  components  in  the  fuselage  by 
adding  higher  harmonic  components  to  the  rotor  controls.  Amplitude  and  phase  of  these 
harmonics  are  determined  on  line  by  a  suitable  control  law.  If  the  higher  harmonic  inputs 
are  applied  in  the  rotating  system,  via  rotating  high  frequency  actuators,  the  technique  is 
usually  called  Individual  Blade  Control  (IBC).  IBC  and  HHC  are  generally  considered  as 
effective  techniques  for  vibration  reduction,  but  issues  of  power  requirements  and  reliability 
have  until  now  prevented  widespread  application  on  production  helicopters  [39]. 

HHC  has  been  the  subject  of  extensive  research  over  the  last  three  decades.  The  research 
until  1982  has  been  reviewed  by  Johnson  [25]  as  part  of  an  extensive  study  of  several  types 
of  HHC  algorithms  and  implementation  techniques  that  remains  relevant  to  this  date.  More 
recent  survey  papers  have  been  written  by  Friedmann  and  Millott  [20]  and  Teves  et  al.  [37]. 
Although  HHC  is  generally  studied  in  the  context  of  rotor  control,  the  basic  HHC  algorithm 
has  also  been  successfully  used  for  fuselage-mounted  active  vibration  reduction  [22], 

Figure  7.1  shows  a  typical  architecture  of  an  HHC  system.  While  this  is  not  the  only 
possible  architecture,  nor  is  it  necessarily  the  best,  it  is  important  for  historical  and  practical 
reasons,  and  has  been  extensively  studied  theoretically  and  experimentally  (e.g.,  [35]).  Some 
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helicopter  outputs,  typically  N/rev  components  of  fuselage  vibrations  are  extracted  through 
a  harmonic  analysis,  fed  to  a  controller  that  computes  appropriate  values  of  the  higher 
harmonic  outputs,  which  are  then  injected  into  the  rotor  controls.  Figure  7.2  depicts  a 
hypothetical  situation  in  which  the  HHC  system  is  inoperative  over  the  first  quarter  of  the 
revolution  of  a  reference  blade,  samples  the  desired  output  over  the  next  quarter  to  build  the 
desired  harmonics  to  attenuate,  computes  the  required  inputs  over  the  next  half  revolution, 
and  updates  the  HHC  inputs  at  the  end  of  the  revolution.  These  steps  are  not  instantaneous, 
but  are  carried  out  over  finite  time  intervals,  and  therefore  cannot  be  rigorously  described 
in  state  space  form. 

This  modeling  problem  is  the  probable  reason  why,  although  the  basic  characteristics  of 
HHC  algorithms  have  been  investigated  extensively  in  the  published  literature,  the  influence 
of  the  aeromechanic  behavior  of  the  entire  helicopter  on  the  performance  and,  especially,  the 
stability  of  the  HHC  system  has  been  typically  ignored.  The  helicopter  dynamics  have  only 
been  considered  indirectly,  through  their  contribution  to  the  gain,  or  T-,  matrix.  Stability 
analyses  have  focused  on  the  stability  of  the  HHC  update  algorithm,  with  or  without  on¬ 
line  identification,  not  on  the  closed  loop  stability  of  the  complete  helicopter  with  the  HHC 
system  turned  on.  Also,  no  study  on  the  potential  interactions  between  an  HHC  system  and 
a  flight  control  system  was  available  in  the  literature.  These  analyses  require  mathematical 
representation  in  state-space  form,  preferably  with  constant  coefficients. 

Some  recent  studies  have  addressed  this  gap.  Cheng  et  al.  [9]  have  extracted  a  linearized, 
state-space,  time-invariant  model  of  the  helicopter,  and  have  used  it  for  a  study  of  the  inter¬ 
action  between  the  HHC  and  the  flight  control  system.  The  harmonic  analyzer  is  essentially 
modeled  as  if  it  was  an  analog  analyzer,  continuously  extracting  the  required  N/rev  vibration 
harmonics.  The  computation  delays  are  modeled  using  Pade  approximants  and  the  HHC 
input  is  updated  continuously.  The  gain  matrix  T  is  fixed.  One  important  conclusion  is  that 
while  the  HHC  has  little  influence  on  the  flight  control  system  and  on  the  handling  qualities, 
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the  reverse  is  not  true.  There  is  a  significant  effect  of  maneuvers  on  vibration,  which  leads 
to  higher  required  crossover  frequencies  for  the  HHC  loop. 

Lovera  et  al.  [28]  have  studied  the  closed  loop  stability  of  a  hingeless  rotor  helicopter 
equipped  with  an  HHC  system.  The  T-matrix  is  constant,  and  the  HHC  algorithm  is  written 
as  a  linear  time- invariant  dynamic  compensator  using  the  technique  developed  by  Wereley 
and  Hall  [40],  extended  to  cover  inputs  and  outputs  at  different  harmonics.  This  compensator 
is  coupled  to  a  high-order,  coupled  rotor- fuselage  model,  and  the  closed  loop  stability  is 
studied  using  Floquet  theory  and  constant  coefficient  approximations.  The  key  results  of 
the  study  are  that:  (i)  the  HHC  system  does  not  degrade  the  aeromechanic  stability,  (ii) 
the  time  constants  of  the  HHC  are  such  that  interactions  with  a  flight  control  system  could 
occur,  and  (iii)  that  the  effects  of  periodic  coefficients  need  to  be  taken  into  account. 

A  gap  in  the  analysis  of  HHC  systems  however  remains  open,  because  the  issues  related 
to  the  actual  discrete-time  implementation  were  not  included  in  Ref.  [28]  and  have  never 
been  studied  in  detail  in  the  literature.  Discrete-time  issues  are  likely  to  play  a  major  role 
in  determining  the  closed  loop  behavior  of  the  system.  In  fact,  the  typical  update  frequency 
for  HHC  control  inputs  is  1/rev,  which  is  comparable  with  the  bandwidth  of  the  dynamics  of 
the  open  loop,  i.e.,  uncontrolled  helicopter  rotor.  Previous  work  on  this  problem  (see,  e.g., 
[27,  28,  40])  tried  to  overcome  this  difficulty  by  developing  continuous-time,  time  domain 
counterparts  of  the  discrete  elements.  On  the  other  hand,  a  more  complete  and  rigorous 
picture  of  the  operation  of  HHC  can  be  obtained  by  looking  at  the  entire  control  loop, 
including  the  coupled  rotor-fuselage  dynamics,  in  discrete  time  rather  than  in  continuous 
time. 

In  the  light  of  the  preceding  discussion,  this  chapter  has  the  following  objectives: 

1.  To  describe  a  typical  discrete-time  HHC  architecture  and  derive  suitable  linearized, 
state-space,  discrete  models  for  all  the  components  present  in  the  control  loop; 
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2.  To  formulate  a  coupled  helicopter-HHC  discrete  model,  and  convert  it  from  a  periodic, 
multiple  sampling  rate  model  to  a  constant  coefficient,  single  sampling  rate  model  using 
time-lifting  techniques;  and 

3.  To  perform  a  closed  loop,  aeromechanic  stability  and  response  analysis  of  the  discrete- 
time,  coupled  helicopter-HHC  model,  and  compare  it  with  the  corresponding  results 
obtained  using  a  continuous-time  model. 

7.2  Helicopter  Model 

The  baseline  simulation  model  used  in  this  study  is  a  nonreal-time,  blade  element  type, 
coupled  rotor-fuselage  simulation  model.  The  model  is  discussed  in  detail  in  Ref.  [38],  and 
only  a  brief  description  will  be  provided  here.  The  fuselage  is  assumed  to  be  rigid  and 
dynamically  coupled  with  the  rotor.  A  total  of  nine  states  describe  fuselage  motion  through 
the  nonlinear  Euler  equations.  Fuselage  and  blades  aerodynamics  are  described  through 
tables  of  aerodynamic  coefficients,  and  no  small  angle  assumption  is  required.  A  coupled 
flap-lag-torsion  elastic  rotor  model  is  used.  Blades  are  modeled  as  Bernoulli-Euler  beams. 
The  rotor  is  discretized  using  finite  elements,  with  a  modal  coordinate  transformation  to 
reduce  the  number  of  degrees  of  freedom.  The  clastic  deflections  are  not  required  to  be 
small.  Blade  element  theory  is  used  to  obtain  the  aerodynamic  characteristics  on  each  blade 
section.  Quasi-steady  aerodynamics  is  used,  with  a  3-state  dynamic  inflow  model.  The  trim 
procedure  is  the  same  as  in  Ref.  [8].  The  rotor  equations  of  motion  are  transformed  into  a 
system  of  nonlinear  algebraic  equations  using  a  Galerkin  method.  The  algebraic  equations 
enforcing  force  and  moment  equilibrium,  the  Euler  kinematic  equations,  the  inflow  equations 
and  the  rotor  equations  are  combined  in  a  single  coupled  system.  The  solution  yields  the 
harmonics  of  a  Fourier  expansion  of  the  rotor  degrees  of  freedom,  the  pitch  control  settings, 
trim  attitudes  and  rates  of  the  entire  helicopter,  and  main  and  tail  rotor  inflow. 
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Linearized  models  are  extracted  numerically,  by  perturbing  rotor,  fuselage,  and  inflow 
states  about  a  trimmed  equilibrium  position.  The  resulting  continuous-time,  linearized, 


time  periodic  model  of  the  helicopter  is  in  the  form 

Xf/(t)  =  +  BFcs(t)llFCs(t)  + 


(7.1) 


y  Hit)  =  Cn(t)-X.H(t )  +  DFcs(t)\lFCs(t )  +  DHHc(t)~U-HHc(t)- 
where  all  the  matrices  are  periodic,  with  period  corresponding  to  N/rev;  the  control  vectors 
u Fcs(t)  and  u uHcif)  are  defined  as  in  Eqs.  (7.4)  and  (7.5),  while  the  outputs  are  the  body 
frame  accelerations  measured  by  the  control  system. 

The  matrices  of  the  linearized  model  are  generated  as  Fourier  series.  For  example,  the 
state  matrix  Ar(xJ>)  is  given  as: 


K 

Ah(‘ 0)  =  Aho  +  ^  {AHkc  cos  kN xp  +  AHks  sin  kNxp)  (7.2) 

k= 1 

with  xp  =  Qt,  and  where  the  matrices  Ar 0,  Ar^,  and  A^ks  are  constant. 

The  control  matrices  Bpcsi'P)  and  Brrq (i/0  are  obtained  in  the  same  Fourier  series  form 
as  Ah (-0),  Eq.  (7.2),  by  assuming  that  the  pitch  control  angle  of  the  i-th  blade  is  given  by 


9i(xp)  =  90  +  9ic  cos  xp  +  6 is  sin  xjj  +  6^c  cos  3 xp  +  93s  sin  3xJj  + 

+9 4c  cos  4:Xp  +  9is  sin  Axp  +  95c  cos  5 x/j  +  95s  sin  5xjj  (7.3) 


where  xpi  —  xp  +  2-ni/N ,  %  —  0, . . .  ,  N  —  1  is  the  azimuth  angle,  and  the  number  of  blades 
IV  =  4.  Therefore,  the  input  harmonics  are  defined  in  the  rotating  system,  but  they  are 
identical  for  each  blade.  This  arrangement  will  be  defined  as  HHC,  although  it  could  also 
fall  under  some  definitions  of  IBC.  Note  that  the  theoretical  development  that  follows  is 
entirely  independent  of  the  specific  configuration  for  actuators  and  sensors. 

The  vector  Uf cs(t)  contains  the  controls  that  would  be  applied  by  the  pilot  or  the  flight 
control  system.  For  the  derivations  of  this  chapter  u Fcs(t)  is  defined  as 

u Fcs(t)  =  [0o(t)  9lc(t)  9ls(t)f  (7.4) 
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The  input  vector  u Fcs(t)  actually  used  in  the  simulations  also  includes  the  tail  rotor  col¬ 
lective  pitch  9t(t).  The  partition  u.HHc(t )  contains  the  harmonics  of  the  HHC  system,  that 

is 

uHHc(t)  —  [#3  c{t)  @3  s(t)  @4  c(t)  94s(t)  95c(t)  95s(t)}T  (7.5) 

The  development  that  follows  only  addresses  the  HHC  control  loop.  The  HHC  analysis 
holds,  under  linearity  assumptions,  regardless  of  whether  the  FCS  loops  are  open  or  closed. 

The  output  vector  y#(t)  can  be  formed  with  any  of  the  three  linear  and  three  angular 
components  of  the  accelerations,  measured  at  one  or  more  points  of  the  airframe  (the  di¬ 
mensions  of  u(t)  and  y (t)  need  not  be  the  same).  The  output  matrices  Cn(t),  DFcs(t )  and 
DHHc(t)  in  Eq.  (7.1)  therefore  depend  on  the  specific  form  of  the  output  vector  y H(t)  (i.e., 
on  the  specific  arrangement  of  sensors),  and  will  be  provided  later  in  the  chapter. 

7.3  Higher  Harmonic  Control 

The  HHC  controller  used  in  the  present  study  is  based  on  a  linear,  steady  state,  frequency  do¬ 
main  representation  of  the  dynamics  of  the  helicopter.  The  vector  u HHc(k)  of  the  harmonics 
of  the  rotor  controls  computed  by  the  HHC  system  is  defined  as 

u mic(k)  =  [93c(k)  93s(k )  9Ac(k)  94s(k)  •••  95c(k )  95s(k)]T  (7.6) 

where  k  is  the  discrete-time  index  associated  with  the  sample  time  at  which  the  control  loop 
operates  (1/rev).  The  vector  y HHcik)  contains  the  N/rev  cosine  and  sine  harmonics  of  the 
fuselage  accelerations,  and  is  defined  as 

y HHc(k)  =  \ylfc(k)  y2Nc{k)  ...  y%c{k )  VNS(k)  yNS(k)  ...  yNs(^)]  (7-7) 
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where  ylNc[k)  and  ylNs(k),  i  —  1, . . .  ,n  are,  respectively,  the  cosine  and  sine  components  of 
the  i-th  N/rev  output  yl(k),  each  defined  as 

i  /‘(k+l)n 

VNc(k)  =  -  VhW  COS  (7.8) 

^  J  kn 

1  l'(k+l)'TT 

VNs(k)  =  -  VhW  sin  N ijdil)  (7.9) 

^  J  k-K 

The  generic  output  yl{k)  can  be  one  of  the  linear  or  angular  components  of  the  accelerations, 
which  in  turn  can  be  measured  at  one  or  more  locations.  In  general,  there  will  be  n  such 
measurements,  and  therefore  the  measurement  vector  will  have  p  =  2 n  elements.  Finally,  y° 
is  the  vector  of  N/rev  disturbance  components  corresponding  to  y HHc(k). 

Then,  assuming  that  the  accelerations  are  linearly  related  to  the  HHC  harmonics,  the 
variables  defined  above  can  be  related  by 

y  HHc(k)  =  TcuHHC(k )  +  y  °(k)  (7.10) 

where  Tc  is  a  real,  constant  coefficient  matrix  that  links  the  harmonics  of  the  HHC  inputs 
to  those  of  the  acceleration  response  (see  Ref.  [28]  for  details).  The  matrix  Tc  can  be 
either  estimated  from  measured  data  using  on  line  or  off  line  identification  algorithms,  or 
computed  on  the  basis  of  a  mathematical  model  of  the  helicopter  as  done  in  the  present 
study.  In  general,  Tc  is  a  function  of  the  flight  condition.  At  each  discrete  time  step  the 
HHC  controller  selects  the  value  of  the  input  harmonics  Uhhc  to  reduce  the  effect  of  y°  on 
y hhc-  Assuming  that  the  baseline  acceleration  vector  y°  is  constant  over  the  time  step,  the 
optimal  open  loop  solution  is  given  by 

u  HHcik)  =  -Tlr(k)  (7.11) 

where  T (t  is  the  pseudoinverse  of  the  Tc  matrix  (which  is  not  necessarily  square).  Since 
in  general  y°  cannot  be  measured  directly,  the  same  result  can  be  obtained  by  using  a 
discrete-time  integral  control  law  in  closed  loop,  i.e.,  based  on  the  measurements  of  y hhc'- 

ilHHc(k  +  1)  =  u HHcik)  —  TcY  HHc{k)  (7-12) 
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The  HHC  control  algorithm  used  in  the  present  study  is  defined  in  terms  of  the  mini¬ 
mization  of  a  quadratic  cost  function  J  of  the  form  [33] 

J  =  y  'hhcW'Yhhc  +  ^hhcR^-hhc  (7-13) 

where  W  =  WT  >  0  and  R  =  R1  >  0  are  matrices  that  allow  different  weighting  of 
acceleration  outputs  and  Au HHc  is  the  increment  of  u HHc  from  one  iteration  to  the  next: 

Au HHc(k  +  1)  =  u  HHc{k  +  1)  —  u HHc(k)  (7-14) 

The  minimization  of  J,  Eq.  (7.13)  leads  to  the  control  law 

A  uHHC(k  +  1)  =  —(TqWTc  +  Ry^Wynncik)  (7.15) 

The  Tc  matrix  links  the  N/rev  harmonics  of  the  output  to  the  harmonics  of  the  HHC  input 
vector  uhhc-  The  matrix  is  fixed,  and  is  obtained  from  the  linearized  model  of  the  complete 
helicopter  using  a  methodology  based  on  the  harmonic  transfer  function.  The  derivation  is 
presented  in  detail  in  Ref.  [28]. 

7.4  Architecture  of  the  HHC  System 

The  present  study  focuses  on  the  simulation  of  the  HHC  architecture  shown  in  Fig.  7.1.  The 
operation  of  the  system  consists  of  the  following  three  steps,  which  are  performed  at  every 
rotor  revolution:  (i)  the  determination  of  the  N/rev  acceleration  output  vector  y HHc(k)',  (ii) 
the  update  UHHc(k  +  l)  of  the  control  input  using  Eq.  (7.15);  and  (iii)  the  actual  application 
of  u HHc(k  +  1)  via  a  simple  zero-order-hold. 

Two  different  sampling  rates  are  used.  The  first,  and  faster,  corresponds  to  the  sampling 
of  the  acceleration  signals  required  to  reconstruct  the  N/rev  harmonics.  The  second,  and 
slower,  is  that  corresponding  to  the  1/rev  update  of  the  HHC  inputs.  The  portions  of  the 
block  diagram  where  each  sampling  rate  is  used  are  shown  in  Fig.  7.1. 

The  closed-loop  analysis  of  the  helicopter  with  the  HHC  system  is  carried  out  as  follows: 
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1.  Discrete-time  models  are  obtained  for  each  of  the  components  of  the  control  loop, 
including  the  block  representing  the  dynamics  of  the  helicopter; 

2.  A  complete  model  is  obtained  for  the  series  connection  of  hold  circuit,  helicopter,  and 
harmonic  analyzer:  this  model  will  prove  to  be  time-periodic; 

3.  A  time-invariant  reformulation  of  the  complete  model  is  obtained  using  the  theory 
of  time-lifting  of  periodic  systems,  using  the  slower  sampling  rate  (i.e.,  that  of  the 
controller) ; 

4.  The  overall  closed-loop  stability  analysis  is  carried  out  in  discrete-time. 

7.5  Discrete  Models  of  the  Loop  Components 

In  this  section,  discrete-time  models  of  all  the  elements  in  the  closed  loop  scheme  of  Fig.  7.1 
are  derived.  They  include:  (i)  the  helicopter  model,  (ii)  the  harmonic  analyzer,  (iii)  the 
controller,  and  (iv)  the  zero-order-hold. 

According  to  the  architecture  defined  in  the  previous  section,  the  HHC  inputs  are  up¬ 
dated  at  1/rev,  while  the  outputs  are  sampled  at  a  higher  frequency  in  order  to  allow  the 
reconstruction  of  the  N/rev  component  of  the  accelerations  of  interest. 

7.5.1  Discrete  helicopter  model 

The  discrete-time  helicopter  dynamic  model  is  obtained  from  the  linearized  continuous-time 
model  of  Eq.  (7.1).  The  sampling  frequency  is  the  faster  of  the  two  in  the  system,  i.e.,  that 
required  to  allow  the  reconstruction  of  the  N/rev  component  of  the  accelerations  of  interest. 

Writing  the  state  and  output  equations,  Eq.  (7.1),  for  the  helicopter  model  over  a  sam¬ 
pling  interval  P,  under  the  usual  assumption  of  constant  input  in  the  interval,  one  can  write 
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the  analytical  solution  for  the  continuous  state  vector  x#  (t)  at  time  t  =  77P  +  P  as 


j-rfP+P 

xH(r)P  +  p)  =  ®h(vP  +  P,  r)P)xH{r]P)  +  /  ®h{vp  +  p,  T)BHHC(T)uHHc(T)dT 

J  r/P 

(7.16) 

y  h(vP)  =  Ch{,tiP)'x-H{vP)  +  DHHC(r]P)uHHC{r]P)  (7-17) 


where  is  the  state  transition  matrix  associated  with  the  state  matrix  Ah-  Defining  the 
discrete-time  state  vector  x#  as  x#  (77)  —  ^h^P)  (similarly  for  the  control  vector  u Hhc  and 
the  output  vector  yH)  results  in  the  following  discrete-time,  state-space  linearized  model  of 
the  helicopter 

xjf (7 7  +  1)  =  AH{rj)x.H{rj)  +  Bhhc(v)uhhc{v)  ^  lgN 

y  h{v)  =  CH(r))5cH{v)  +  DHHC(r))uHHc(ri ) 

where  the  system  matrices  are  defined  as 


Ah(v)  =  ®h(vP  +  P,VP) 

Bh(v )  =  f  ®H(rjP  + P,t' +  rjP)BHHc(r' +  rjP)dr' 
Jo 

CH(V)  =  CH(VP) 

Dhhc{i ?)  =  DHHC(rjP ) 


(7-19) 

(7.20) 

(7.21) 

(7.22) 


The  system  matrices  are  periodic,  with  a  common  period  equal  to  ns  samples. 

7.5.2  Harmonic  Analysis 

To  implement  the  HHC  control  algorithm,  the  N/rev  components  y hhc  of  the  output  accel¬ 
erations  must  be  extracted  from  their  time  domain  measurements  y H-  hi  each  period,  the 
information  about  y#  is  available  starting  from  77  =  nsj 2  —  1  +  Kns,  K  =  1,2, .. .  but  is 
provided  as  output  only  at  77  =  (. K  +  1  )ns.  The  operation  of  the  harmonic  analyzer  can  be 
described  mathematically  by  a  linear  time-periodic  model  with  discrete  time  77  and  period 
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xF(r 7  +  1)  =  AF(rj)5cF(rj)  +  BF{rj)yH{rj) 


(7.23) 


ns: 


y  HHcij])  =  CF(rj)yLF{rj) 

The  matrices  AF(p),  BF(r /)  and  CF(g)  are  defined  in  a  way  that  reflects  the  various  phases 
of  the  harmonic  analysis  that  occur  over  one  rotor  revolution.  For  the  four-bladed  rotor  of 
this  study  there  are  three  distinct  phases,  defined  as  follows  (see  also  Figure  7.2): 

1.  During  the  first  quarter  of  the  period ,  rj  =  1,2,...  ,ns/ 4,  the  output  signal  yF  is 
allowed  to  reach  steady  state  following  the  update  of  the  control  input  u hhc  at  the 
end  of  the  previous  rotor  revolution.  During  this  time,  the  output  of  the  harmonic 
analyzer  is  set  to  zero  and  the  y#s  measured  are  not  accumulated  in  the  integrals, 
Eqs.  (7.8)  and  (7.9).  Therefore  the  state  space  matrices  of  the  harmonic  analyzer  are 
given  by 

Mv)  =  I  BF(V)  =  0  CF(ri)  =  0  (7.24) 


2.  During  the  second  quarter  of  the  period ,  rj  =  ns/4+l,  nsf  4+2, . . .  ,  ns/ 2,  the  outputs  yF 
are  actually  sampled  and  the  integrals,  Eqs.  (7.8)  and  (7.9),  computed,  but  the  output 
of  the  analyzer  is  still  kept  to  zero.  Note  that  the  state  of  the  harmonic  analyzer  must 
be  reset  to  zero  at  the  beginning  of  the  second  quarter  of  the  period,  so  the  state  space 
matrices  will  be  defined  as: 

Mri)  =  {  °  ’’ =  T  (7.25) 

1  elsewhere 


Bf(v) 


blkdiag  < 


(  2Ntt 
sin  - rj 

V  n8 


2Nn 

cos  |  - rj 

ns 


Cf(v )  = 


0 
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(7.26) 

(7.27) 


3.  During  the  remaining  half  period ,  77  =  ns/2  +  1,  ns/2  +  2, . . .  ,ns  the  new  value  of  the 
control  input  u HHC  is  computed  and  applied  to  the  rotor;  therefore  the  operation  of 
the  analyzer  is  stopped  and  the  computed  value  for  the  N/rev  harmonic  of  y#  is  made 
available  at  the  end  of  the  period: 


Mv)  =  1 

(7.28) 

BF{rj)  =  0 

(7.29) 

r  8  r 

—I  rj  =  ns 

cF(v)  =  Us 

(7.30) 

(  0  elsewhere 

The  entire  sequence  of  operations  described  above  can  be  summarized  in  the  following 
expressions  for  the  state  space  matrices  of  the  harmonic  analyzer,  which  hold  for  a  generic 
value  of  the  1/rev  discrete  time  index  k  (i.e. ,  for  the  generic,  /c-th  rotor  revolution): 


Mv)  = 


Tic 


0  rj  —  kns  +  — 
/  elsewhere 


(7.31) 


ns  ns 

Mv)  =  i  1  ‘"•  +  Ts,<b,-+2 


0 


elsewhere 


(7.32) 


Bf(v)  =  PF(rj)b1kdi&g  < 


2Nn 

sm  |  - rj 

ns 


cos 


CF(v)  = 


— /  r)  =  kns  +  ns 
ns 

0  elsewhere 


(7.33) 

(7.34) 


Therefore,  the  output  of  the  above  model  is  nonzero  only  for  77  =  0,  ns,  2 ns, . 
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7.5.3  Controller 


The  control  law  given  by  Eq.  (7.15)  can  be  written  in  state-space  form  as  a  linear  time- 
invariant  system 


xc(fc  +  1)  =  Acxc(k)  +  BcyHHcik)  (7.35) 

u  HHc(k)  =  CcZc(k)  (7.36) 

where 

Ac  =  I  (7.37) 

Bc  =  -{t£WTc  +  R}-1T£w  (7.38) 

Cc  =  I  (7.39) 

7.5.4  Zero-Order-Hold  circuit 


The  hold  circuit  is  the  interface  between  the  controller  and  the  helicopter.  Since  the  controller 
operates  at  the  discrete-time  k  (i.e.,  once  per  revolution)  while  the  helicopter  model  has  been 
obtained  at  the  discrete-time  rj  (i.e.,  once  per  sample  needed  to  extract  the  N/rev  harmonics), 
the  controller  provides  a  new  value  of  the  control  variables  only  at  r)  =  kns ,  k  =  1,2, ... , 
and  this  output  must  be  kept  constant  for  the  intervening  samples  kns  <  rj  <  (k  +  l)ns. 
Therefore,  the  model  of  the  hold  circuit  is  linear,  discrete-time  periodic,  with  discrete  time 
r)  and  period  ns ,  and  is  given  by 


xz(t7  +  1)  =  Az{ri)^z{ji)  +  Bz(rj)uHHc{r))  (7.40) 

u hhc{v)  =  Cz(r])±z(r])  +  Dz(rj)uHHc{ri)  (7-41) 
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where 


and 


Az(  V) 

=  S(v) 

(7.42) 

Bz(n) 

=  I-S(v) 

(7.43) 

Czivl 

=  S(rj) 

(7.44) 

bz(n) 

=  I~S(r}) 

(7.45) 

s(v) 


0  T)  —  kns ,  k  —  1,2,... 
/  elsewhere 


(7.46) 


7.5.5  Series  connection  of  zero  order  hold,  helicopter  model  and  harmonic  an¬ 
alyzer 


The  overall  discrete  HHC  model,  which  relates  the  harmonics  of  the  HHC  input  u hhc  to 
the  harmonics  of  the  acceleration  output  y hhc  can  be  obtained  by  connecting  in  series  the 
harmonic  analyzer,  Eq.  (7.23),  the  discrete  model  for  the  response  of  the  helicopter  to  HHC 
inputs,  Eq.  (7.16),  and  that  of  the  zero  order  hold,  Eq.  (7.40).  The  model,  with  discrete-time 
r),  is  given  by 


xz(r?  +  1) 

=  Azx-z(r])  +  BzxiHHc(;ri) 

(7.47) 

xjf  (r 7  +  1) 

=  AH5tH(;ri)  +  BHHcCz*z(:n)  +  BhhcDzUhhc(v) 

(7.48) 

xf() 7  +  1) 

=  AF5tF(r])  +  BfCh±h{  v)  +  BFbHHCCz^-z{ri)  + 

+BFbHHcbzUHHc(v) 

(7.49) 

y  HHciv) 

=  CFiF(rj) 

(7.50) 

(the  argument  r]  in  the  matrices  has  been  omitted  for  simplicity).  This  model  cannot  be 
connected  directly  to  the  HHC  controller  because  its  sampling  rate  is  still  different  from  that 
of  the  discrete  HHC  control  law  (ns/ rev  vs.  1/rev).  A  combined  model  at  the  same  sampling 
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rate  as  the  HHC  controller  can  be  obtained  using  the  theory  of  time-invariant  reformulations 
of  linear  time-periodic  systems  and,  more  precisely,  through  a  time-lifted  reformulation. 

7.6  Time-lifted  Formulation  of  the  Closed  Loop  System 

A  closed-loop  stability  analysis  requires  a  single  sampling  period  for  the  entire  system.  There¬ 
fore,  the  model  will  be  reformulated  using  the  shorter  sampling  period  P  =  T /ns  as  the 
common  period.  The  implicit  assumption  that  the  longer  sampling  period  T,  correspond¬ 
ing  to  one  rotor  revolution,  is  an  integer  multiple  of  the  acceleration  sampling  period  P  is 
clearly  quite  reasonable.  This  reformulation  is  carried  out  using  time  lifting,  and  results  in 
a  discrete  model  with  overall  sampling  period  T  and  time  index  k.  Time  lifting  is  described 
in  detail  in  Ref.  [4]  and  is  summarized  in  Ref.  [11],  which  also  contains  numerical  examples 
of  the  application  to  a  flapping  rigid  rotor  blade.  A  useful  byproduct  of  the  use  of  a  lifted 
reformulation  is  that  the  resulting  model  has  constant  coefficients,  which  simplifies  the  closed 
loop  stability  analysis. 

7.6.1  Time  lifting  of  periodic  systems 

Time  lifting  is  based  on  the  idea  that  the  knowledge  of  the  state  vector  at  time  k  and  of 
the  inputs  between  time  k  and  k  +  1  is  sufficient  to  determine  the  value  of  the  state  at  time 
k  T  1  and  the  value  of  the  outputs  between  k  and  k  +  1.  The  key  steps  of  time  lifting  will 
be  briefly  outlined  here. 

Consider  the  linear,  discrete-time  periodic  system  with  a  period  equal  to  ns  samples,  and 
with  m  inputs  and  p  outputs 

x(f+l)  =  A(t)k(t)  +  B(t)u{t)  .  , 

y  (t)  =  C(t)k(t)  +  D(t)u(t) 

where  t  =  0,1,2,...  ,ns  —  1.  Then,  the  state  vector  at  time  t  >  t  is  given  by  the  discrete-time 
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Lagrange  formula  [34] 


x(i)  =  'L(£,t)x(t)  +  ~  !)u(j  -  !)• 

j=r+l 


(7.52) 


where  \k(£,  r)  =  7l(f  —  l)A(t  —  2)  •  •  •  A(t)  is  the  transition  matrix  from  time  r  (also  an  integer 
between  0  and  ns  —  1)  to  time  t  for  the  state  equation,  Eq.  (7.51).  Equation  (7.52)  can  be 
used  to  build  an  equivalent  (i.e.,  with  the  same  output  given  the  same  input)  time-invariant 
system  by  sampling  the  state  vector  at  a  frequency  of  ns ,  and  packing  the  input  and  output 
vectors  for  each  sample  into  larger  input  and  output  vectors.  This  results  in  the  “lifted” 
reformulation  [4,  11] 


x(/c  +  1)  =  Fx(fe)  +  Guuft(k) 
y  uft(k)  =  HZ(k)  +  EuHft(k) 

where  the  extended  input  vector  uuft(k )  has  size  mns  and  is  defined  as 

u uftik)  =  [u {kns)T  ■  ■  ■  u (kns  +  ns  -  l)T]T 
and  the  extended  output  vector  y uft{k)  has  size  pns  and  is  defined  as 


(7.53) 


y uft(k)  =  [  y (kns 


y(kns  +  ns  —  1)T 


The  time  invariant  system  matrices  F,  G,  H,  and  E  have  dimensions,  respectively,  of  n 
by  n,  n  by  mns,  pns  by  n,  and  pns  by  mns  and  are  given  by  [4,  11] 


F  = 

A(ns  —  l)A(ns  —  2)  •  •  •  A(0) 

(7.54) 

G  = 

[T(ns,l)fi(0)  T(ns,  2)5(1)  'ff(ns,ns)B(ns  -  1)] 

(7.55) 

H  = 

[C(0f  T(l,0)TC(lf  T(ns-l,0)rCK-lf]T 

(7.56) 

E  = 

L  j  =  I;  2, . . .  ns 

(7.57) 

with 

(  0  %  <  j 

Eij  =  <  D(i  -  1)  i  —  j 

[  C(i  -  1 W  -  1  -  1)  i  >  j 

(7.58) 
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It  should  be  noted  that  the  matrix  F  is  the  Floquet  Transition  matrix  of  the  discrete-time 
periodic  system,  Eq.  (7.51),  therefore  the  constant-coefficient  lifted  system  of  Eq.  (7.53) 
has  exactly  the  same  stability  characteristics  as  the  time-periodic  system.  This  Floquet 
Transition  Matrix  is  also  the  same  as  that  of  the  original  continuous  system,  and  therefore 
continuous  and  discrete  systems  also  have  the  same  stability  characteristics  [4,  11], 

7.6.2  Lifted  form  of  the  HHC  loop 

Time  lifting  can  be  directly  applied  to  the  coupled  hclicopter/HHC  system,  Eqs.  (7.47)- 
(7.50).  The  corresponding  lifted  form  is  given  by: 


x(/c  +  1)  =  Fx(fc)  +  Guuft(k) 

y  uftik)  =  H±(k)  +  Euuftik) 


The  input  vector  u uft(k)  is  given  by 

uuft(k) 


I 

I 


iiHHc(k) 


I 


(7.59) 

(7.60) 


(7.61) 


where  the  mns  by  m  matrix  that  relates  u uft(k)  and  u HHc(k)  is  obtained  by  assembling  ns 
identity  matrices  because  the  HHC  input  vector  u HHcik)  is  held  constant  over  all  the  ns 
samples  that  make  up  one  rotor  revolution. 

Similarly,  once  the  lifted  output  vector  yuft  has  been  obtained  as  the  output  of  the  lifted 
system,  Eqs.  (7.59)-(7.60),  the  actual  discrete  output  vector  y HHcik)  can  be  recovered  by 
observing  that  it  is  the  output  of  the  Fourier  coefficients  extractor,  which  is  only  evaluated 
once  per  revolution,  i.e.,  at  times  rj  =  0 ,ns,  2 ns, ....  Therefore 


y HHc(k)  =  [I  0]y uft(k)  (7.62) 

where  y HHc(k)  has  size  p  and  y uft(k)  has  size  pns.  On  the  basis  of  Eqs.  (7.61)  and  (7.62), 
and  recalling  the  state  space  form  for  the  discrete  HHC  controller,  Eqs.  (7.35)-(7.36)  the 
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overall  closed  loop  system  can  be  constructed  as  follows: 


x(/c  +  1) 


xc(&  +  1) 


I 

I 


F5t(k)  +  G 


Cc*c(k) 


Bc  [I  0 ]H±{k)  + 


/ 

'  I ' 

\ 

I 

AC  +  [I  0 }E 

V 

I 

) 

*c{k) 


(7.63) 


(7.64) 


This  closed  loop  helicopter/HHC  system  is  now  linear  time-invariant,  with  discrete-time 
k  and  state  variables  xz,  xh,  %f,  and  xc ■  Therefore,  the  closed  loop  stability  analysis  can 
be  carried  out  by  checking  the  eigenvalues  of  the  closed  loop  system. 


7.7  Results 


This  section  presents  closed-loop  stability  and  response  results  for  a  coupled  helicopter- 
HHC  system.  The  stability  results  and  the  linearized  time  histories  are  obtained  from  the 
linearized,  time-lifted  model  of  Eqs.  (7.63)  and  (7.64).  The  closed-loop  response  results  are 
obtained  from  the  full  nonlinear  simulation  of  the  coupled  helicopter-HHC  system. 

The  helicopter  configuration  used  for  the  present  study  is  similar  to  the  Eurocopter  B0- 
105,  with  a  thrust  coefficient  Ct/ct  =  0.071.  Three  blade  modes  are  used  in  the  modal 
coordinate  transformation,  namely,  the  fundamental  flap,  lag,  and  torsion  modes,  with  nat¬ 
ural  frequencies  of  1.12/rev,  0.7/rev,  and  3.4/rev,  respectively.  Because  the  aerodynamic 
model  consists  of  a  simple  linear  inflow  with  quasi-steady  aerodynamics,  vibratory  loads  and 
CG  accelerations,  and  consequently  also  HHC  inputs,  tend  to  be  underestimated.  Therefore, 
their  absolute  values  can  be  considered  representative  only  in  a  qualitative  sense.  However, 
the  overall  simulation  model  is  likely  reasonable  for  stability  studies,  and  for  a  general  as¬ 
sessment  of  the  design  and  closed-loop  analysis  methodology. 

For  all  the  vibratory  response  results,  the  helicopter  is  first  trimmed  in  steady,  straight 
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flight  at  the  desired  velocity  with  the  HHC  system  turned  off.  Then,  the  nonlinear  simulation 
begins,  with  the  pilot  controls  held  fixed  at  their  trim  values  and  the  HHC  system  turned 
on  at  time  t  —  0. 

The  controller  has  been  implemented  in  discrete  time  using  a  ’’fast”  sampling  rate  of 
16  samples  per  rotor  revolution.  The  weighting  matrices  W  and  R  which  define  the  HHC 
quadratic  performance  index,  Eq.  (7.13),  have  been  chosen  to  be  proportional  to  the  identity 
matrix,  i.e.,  W  =  wl  and  R  =  rl. 

7.7.1  Results  for  V=80  kts 

Figure  7.3  shows  the  closed-loop,  vertical  acceleration  response  w  for  three  different  values 
of  the  tuning  parameter  r,  namely  r  —  2  •  104  (top  plot),  r  —  5  •  104  (center  plot),  and 
r  =  105  (bottom  plot),  corresponding  to  increasing  restrictions  on  the  control  effort  or, 
equivalently,  to  decreasing  gain.  Note  that  the  scales  on  the  vertical  axis  are  different  for 
each  plot.  For  the  value  r  —  2  •  104  the  response  loosely  resembles  a  limit  cycle,  but  the 
values  of  the  accelerations  are  of  up  to  ±1  g,  and  therefore  unacceptably  high.  For  this  value 
of  r  the  linearized  analysis  predicts  an  instability.  For  the  value  r  —  5  •  104  the  response 
is  unstable,  but  the  values  of  w  are  now  considerably  smaller,  and  remain  within  ±0.04g 
within  the  first  seven  seconds.  The  linearized  analysis  indicates  a  mild  instability.  Finally, 
for  r  =  105,  the  response  becomes  stable  and  the  HHC  is  very  effective  in  suppressing  w. 
This  is  also  predicted  by  the  linearized  analysis.  The  behavior  of  roll  and  pitch  accelerations 
is  qualitatively  very  similar  to  that  of  w,  and  it  will  not  be  shown  here. 

Similar  plots  had  been  presented  in  Ref.  [28]  for  the  same  flight  condition  and  helicopter 
configuration,  but  with  the  HHC  closed  loop  modeled  as  entirely  continuous,  and  the  har¬ 
monic  analysis,  the  computation  delay,  and  the  zero-order-hold  not  modeled  at  all.  For 
values  of  the  tuning  parameter  r  ranging  from  0  to  104  the  response  was  stable,  and  the 
HHC  would  effectively  suppress  vibrations  in  10  seconds  or  less.  Comparing  these  results 
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with  those  shown  in  Fig.  7.3,  it  is  clear  that  the  allowable  gains  of  the  HHC  system  must 
be  lower  because  of  the  reduction  in  phase  margin  brought  about  by  the  delays  in  the  HHC 
loop. 

Selected  closed  loop  stability  eigenvalues  are  shown  in  Fig.  7.4  as  a  function  of  r.  In 
the  top  portion  of  the  figure  the  eigenvalues  are  presented  in  root  locus  plot  form,  in  the 
bottom  portion  the  real  parts  of  the  eigenvalues  are  plotted  as  a  function  of  r.  These  are  the 
eigenvalues  of  the  discrete  system,  converted  to  continuous  form.  The  closed  loop  system 
becomes  unstable  for  r  <  1.5  •  105.  This  instability  was  not  captured  by  the  approximate 
continuous  HHC  model  of  Ref.  [28],  which  instead  predicted  closed  loop  stability  for  every 
value  of  r.  This  clearly  demonstrates  that  neglecting  the  discrete  nature  of  the  HHC  loop  is 
unconservative,  and  should  be  avoided. 

Figure  7.5  shows  the  time  history  of  the  response  of  just  the  4/rev  harmonic  of  w. 
The  three  curves  show,  respectively,  the  baseline  response  with  the  HHC  system  turned 
off,  the  response  predicted  by  the  nonlinear  simulation  model,  and  that  predicted  by  the 
linearized,  time  periodic  model  used  to  design  the  HHC  system.  Apart  from  a  small  initial 
transient,  caused  by  a  small  initial  mismatch  between  the  trimmed  and  the  time- marching 
solution,  the  baseline  4/rev  response  rapidly  converges  to  a  constant,  nonzero  steady  value. 
The  nonlinear  closed-loop  response  exhibits  a  brief  but  strong  transient,  during  which  the 
acceleration  increases  by  almost  three  times  the  baseline  value.  The  transient  lasts  for  less 
than  2  seconds,  after  which  the  4/rev  response  is  rapidly  reduced  to  almost  zero.  This 
strong  transient  is  not  present  in  the  p  and  q  4/rev  responses,  not  shown  in  the  figure,  which 
start  being  attenuated  as  soon  as  the  HHC  is  turned  on.  The  figure  also  shows  that  the 
4/rev  responses  predicted  by  the  LTP  and  by  the  nonlinear  model  are  nearly  identical.  This 
indicates  that,  for  the  type  of  mathematical  model  used  in  this  study,  and  for  the  flight 
condition  considered,  (i)  the  effect  of  nonlinearities  on  the  4/rev  w  response  is  small,  and  the 
response  is  adequately  captured  by  a  linearized  model  as  long  as  the  periodicity  is  retained 
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(the  same  conclusions  hold  for  p  and  q),  and  (ii)  the  LTP  model  is  sufficiently  accurate  for 
the  HHC  design. 

The  magnitude  of  the  3/,  4/,  and  5/rev  control  harmonics  are  plotted  in  Fig.  7.6  as  a 
function  of  time.  The  controls  for  both  the  rigorous  discrete  model  and  the  approximate 
continuous  model  are  shown  in  the  figure.  Discrete  and  continuous  controls  tend  to  the  same 
magnitude  value,  but  the  initial  transients  are  quite  different.  The  magnitudes  of  the  discrete 
controls  grow  much  faster  than  those  of  the  continuous  controls.  The  magnitudes  of  the  3/ 
and  4/rev  harmonics  of  the  discrete  controller  reach  essentially  the  respective  steady  state 
values  in  about  2  seconds,  that  of  the  5/rev  harmonic  in  about  6  seconds.  The  respective 
values  for  the  continuous  case  are  well  over  6  seconds  for  the  3/  and  5/rev,  and  about  6 
seconds  for  the  4/rev.  This  behavior,  which  is  to  some  extent  counterintuitive  given  that 
the  discrete  HHC  is  operating  at  a  lower  gain  than  the  continuous  HHC,  is  responsible  for 
the  faster  vibration  reduction  of  the  discrete  HHC. 

Better  agreement  between  the  phases  of  the  discrete  and  continuous  model  can  be  seen 
in  Fig.  7.7.  Except  for  the  first  1-2  rotor  revolutions  (slightly  more  for  the  4/rev  control), 
discrete  and  continuous  phases  are  almost  identical. 

7.7.2  Results  for  V=140  kts 

Figure  7.8  shows  the  closed-loop,  vertical  acceleration  response  w  for  three  different  values 
of  the  tuning  parameter  r,  namely  r  =  104  (top  plot),  r  —  5  •  104  (center  plot),  and  r  =  105 
(bottom  plot).  The  scale  on  the  vertical  axis  of  the  top  plot  is  different  from  those  of 
the  other  two.  As  in  the  80  kts  case,  for  the  value  r  =  104  the  acceleration  response  is 
very  high  and  erratic,  with  peaks  well  over  1  g  in  absolute  value,  and  loosely  resembling  a 
limit  cycle.  For  this  value  of  r  the  linearized  analysis  predicts  an  instability.  For  the  value 
r  —  5  •  104  the  response  is  stable,  and  slowly  decreases  in  magnitude.  Finally,  for  r  =  105, 
the  response  is  stable  and  the  vibrations  are  reduced  very  quickly,  in  less  than  2  seconds 
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from  the  application  of  HHC.  The  linearized  analysis  also  indicates  that  these  last  two  cases 
are  stable.  The  behavior  of  roll  and  pitch  accelerations  is  qualitatively  very  similar  to  that 
of  w,  and  it  will  not  be  shown  here. 

Compared  with  the  corresponding  plots  of  Ref.  [28],  for  the  same  flight  condition  but 
with  a  continuous  HHC  model,  the  allowable  gains  are  lower.  In  Ref.  [28]  the  closed-loop 
system  was  studied  for  values  of  r  ranging  from  0  to  1000,  and  in  all  cases  it  was  stable. 

Figure  7.9  is  similar  to  Fig.  7.5,  but  refers  to  V  —  140  kts.  As  in  Fig.  7.5,  after  a  small 
initial  transient,  the  baseline  4/rev  response  rapidly  converges  to  a  constant,  nonzero  steady 
value.  The  nonlinear  closed-loop  response  no  longer  exhibits  the  strong  transient  observed 
at  V  =  80  kts,  and  the  4/rev  response  is  substantially  reduced  after  just  one  second.  The 
same  happens  for  the  p  and  q  4/rev  responses,  not  shown  in  the  figure.  As  in  the  V  =  80 
kts  case,  the  4/rev  responses  predicted  by  the  LTP  and  by  the  nonlinear  model  are  nearly 
identical,  and  the  same  conclusions  on  the  effect  of  nonlinearities  and  the  adequacy  of  the 
LTP  model  for  design  purposes  apply. 

The  magnitude  of  the  3/,  4/,  and  5/rev  harmonics  are  plotted  in  Fig.  7.10  as  a  function 
of  time.  The  controls  for  both  the  discrete  and  the  continuous  model  are  shown  in  the  figure. 
As  in  the  V  =  80  kts  case,  discrete  and  continuous  controls  tend  to  the  same  magnitude 
value,  although  the  initial  transients  are  quite  different  and  the  magnitudes  of  the  discrete 
controls  grow  much  faster  than  those  of  the  continuous  controls.  The  magnitudes  of  all 
the  harmonics  of  the  discrete  controller  reach  their  respective  steady  state  values  in  about 
2  seconds.  The  respective  values  for  the  continuous  case  are  of  6  seconds  and  more.  An 
almost  perfect  agreement  between  the  phases  of  the  discrete  and  continuous  model  can  be 
seen  in  Fig.  7.11.  Except  for  the  first  few  rotor  revolutions,  discrete  and  continuous  phases 
are  almost  identical. 
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7.7.3  Other  considerations 


The  results  presented  in  this  section  underscore  the  importance  of  a  correct  modeling  of 
“real  life”  effects  such  as  discrete  sampling,  computations,  and  control  updates,  even  for 
the  simplified,  fixed  T-matrix  scheme  used  in  this  study.  Several  additional  effects  were 
neglected,  and  should  be  included  or  more  carefully  analyzed  in  future  research.  First,  in  the 
scheme  of  Fig.  7.2,  it  was  assumed  that  the  transient  following  each  HHC  update  would  take 
one  quarter  of  a  rotor  revolution  to  die  out.  Simulation  results  not  presented  here  indicate 
that  a  more  realistic  figure  is  1-2  rotor  revolutions  for  well-damped  rotors  with  mechanical 
lag  dampers,  and  up  to  4-6  revolutions  for  lowly  damped  hingeless  rotors.  Second,  the  HHC 
update  at  each  rotor  revolution  was  simulated  as  a  pure  step.  While  rotating  system  HHC 
actuators  are  very  fast,  they  cannot  generate  such  steps,  and  therefore  they  add  their  own 
delay.  Third,  perfect  measurements  were  assumed,  whereas  real  sensors  introduce  their  own 
dynamics  in  the  loop.  Finally,  practical  digital  harmonic  analyses  will  require  the  use  of 
windows,  which  may  introduce  further  delays  and  spurious  dynamics.  All  these  effects  must 
be  taken  into  account  to  obtain  realistic  predictions  of  the  maximum  performance  achievable 
by  an  HHC  system. 

7.8  Conclusions 

This  chapter  presented  an  aeromechanical  closed  loop  stability  and  response  analysis  of  a 
hingeless  rotor  helicopter  with  an  HHC  system  for  vibration  reduction.  The  analysis  fully 
included  the  rigid  body  dynamics  of  the  helicopter  and  the  flexibility  of  the  rotor  blades.  It 
was  assumed  that  the  gain  matrix  T  was  fixed  and  computed  off-line.  The  discrete  elements 
of  the  HHC  control  loop  were  rigorously  modeled,  including  the  presence  of  two  different 
time  scales  in  the  loop.  By  also  formulating  the  coupled  rotor-fuselage  dynamics  in  discrete 
form,  the  entire  coupled  hclicopter-HHC  system  could  be  rigorously  modeled  as  a  discrete 
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system.  Finally,  the  effect  of  the  periodicity  of  the  equations  of  motion  was  rigorously  taken 
into  account  by  converting  the  system  with  periodic  coefficients  into  an  equivalent  system 
with  constant  coefficients  and  identical  stability  properties  using  a  time  lifting  technique. 

The  most  important  conclusion  of  the  present  study  is  that  the  discrete  elements  in 
the  HHC  loop  must  be  modeled  in  any  HHC  analysis.  Not  doing  so  is  unconservative. 
For  the  helicopter  configuration  and  HHC  structure  used  in  this  study,  an  approximate 
continuous  modeling  of  the  HHC  system  indicated  that  the  closed  loop,  coupled  hclicopter- 
HHC  system  was  always  stable,  whereas  the  more  rigorous  discrete  analysis  shows  that  closed 
loop  instabilities  can  occur.  The  HHC  gains  must  be  reduced  to  account  for  the  loss  of  gain 
margin  brought  about  by  the  discrete  elements. 

Other  conclusions  of  the  study  are:  (i)  the  HHC  is  effective  in  quickly  reducing  vibrations, 
at  least  at  its  design  condition;  (ii)  a  linearized  model  of  helicopter  dynamics  is  adequate 
for  HHC  design,  as  long  as  the  periodicity  of  the  system  is  correctly  taken  into  account,  i.e., 
periodicity  is  more  important  than  nonlinearity,  at  least  for  the  mathematical  model  used 
in  this  study;  and  (iii)  when  discrete  and  continuous  systems  are  both  stable,  the  predicted 
HHC  control  harmonics  are  in  good  agreement,  although  the  initial  transient  behavior  can 
be  considerably  different. 


156 


Figure  7.1:  Block  diagram  of  a  HHC  loop  (Thin  lines:  continuous  signals;  thick  lines:  ’’fast 
sampling”  discrete  signals;  double  lines:  1/rev  sampling  signals.) 
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Figure  7.2:  Operation  of  the  control  system  over  one  rotor  revolution  period  T,  as  a  function 
of  the  azimuth  angle  0. 
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Figure  7.4:  Selected  closed  loop  stability  eigenvalues  for  V=80  kts  (/i  =  0.188)  as  a  function 
of  the  controller  tuning  parameter  r;  root  locus  plot  (top)  and  real  parts  only  (bottom). 
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Figure  7.5:  Closed  loop  4/rev  vertical  acceleration  response  w  at  the  helicopter  center  of 
mass  for  V=80  kts  (/i  =  0.188);  baseline  open-loop  response,  and  prediction  with  linear  and 
nonlinear  simulation  model. 
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Figure  7.8:  Vertical  accelerations  w  at  the  helicopter  center  of  mass  for  V=140  kts  (/i 
0.330)  and  tuning  parameter  r  =  104  (top),  r  —  5  •  104  (center),  and  r  =  105  (bottom). 
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Figure  7.9:  Closed  loop  4/rev  acceleration  response  w  at  the  helicopter  center  of  mass  for 
V=140  kts  {/i  =  0.330);  baseline  open- loop  response,  and  prediction  with  linear  and  nonlinear 
simulation  model. 
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Figure  7.10:  HHC  control  input  magnitude  in  degrees  for  continuous  and  discrete  models, 
V=140  kts  (n  =  0.330);  3/rev  (top),  4/rev  (center),  5/rev  (bottom). 


Figure  7.11:  HHC  control  input  phase  in  degrees  for  continuous  and  discrete  models,  V=140 
kts  {jjL  =  0.330);  3/rev  (top),  4/rev  (center),  5/rev  (bottom). 
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